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ABSTRACT 


A quantitative analysis was carried out to determine the 
stresses present in the leading-edge segment of a P-3C 
aircraft operating within and outside the normal operating 
envelope of the aircraft. The purpose of the analysis was to 
ascertain whether a specific weakness may exist in the 
leading-edge structure which might endanger future operating 
flight crews. A three-step process consisting of a static 
aeroelastic span-load analysis, an inviscid two-dimensional 
panel method, and finite-element analysis was employed in the 
course of the evaluation. Lift-coefficient distributions from 
the wing span-load analyses were used in the two-dimensional 
panel method to determine the pressure distribution around the 
leading edge, which was then used as input to the finite- 
element analysis. Additionally, static aeroelastic-derived 
wing-twist effects were included in the structural model. The 
results of the analysis suggest that the leading edge segment 
Studied may experience stress levels sufficient to cause 


failure within the normal operating envelope. 
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I. INTRODUCTION 


On 13 February 1988, aU.S. Navy (USN) P-3C Update I Orion 
aircraft experienced an in-flight failure of a wing leading- 
edge section during a high-speed, low-altitude maneuver. The 
aircraft was able (with some difficulty) to return to its 
departure point and make a safe landing. Some three years 
later, on 27 April 1991, an Orion operated by the Royal 
Australian Air Force (RAAF) suffered a similar but more 
extensive in-flight failure when it lost three wing leading- 
edge sections during a maneuver similar to the U.S. Navy 
mishap. The Australian crew was unable to return to the 
runway due to the loss of lift available from the damaged 
wings. The aircraft impacted the ocean surface short of its 
island destination in a nose high-attitude with maximum power 
set on all four engines. One crewmember was killed when a 
propeller tore loose from its engine and entered the fuselage 
where he was sitting. 

The leading edge is the rounded front’ portion of the wing 
which is physically attached to the front wing spar and is 
essential in the production of lift by an airfoil. Loss of a 
.leading-edge segment results in a dramatic reduction in wing 


lift capability combined with a corresponding increase in 


drag. The P-3 has three of these leading-edge segments on 
each wing, separated by the engine nacelles. They are best 
referred to as inboard, center and outboard sections. Figure 
1 shows the P-3 aircraft. The inboard sections are those 
between the fuselage and the 
inboard engines on this 
four-engine Lr Ui ay OO Oya GLSy ©: 
aircraft. The center 
sections (blackened) span 


the distance between the 





inboard and outboard engine | 
nacelles, while the outboard Piguresl. P+3-leadingagia 
sections cover the remaining distance from the outboard 
engines to the wing tips. The length (fore and aft) of these 
leading-edge sections is 15% of the chord (total fore and aft 
distance) of the wings. 

The USN P-3 lost the starboard wing center leading-edge 
section while the RAAF aircraft lost both of its center 
sections and the starboard inboard section. This study 
focuses on the center leading edge sections. 

While both of these mishaps. were investigated by the 
appropriate authorities, the cause of the failures were 
undetermined, though widely suspected to be the result of 
aircraft overstress due to pilot error. In both cases, the 


pilots did not feel that they had exceeded the normal 


operating envelope (this envelope will be detailed later) for 


the aircraft. Even if the aircraft were operated outside the 
authorized envelope, the question remains as to why the 
leading edge sections failed before some other component. 
This thesis was undertaken for the purpose of studying the 
aerodynamic loading and structural response of a leading edge 
section due to operation within and outside the normal flight 
envelope. Ultimately, its aim was to determine whether there 
may be a need to further restrict the operating envelope or 
recommend some modification to the existing structure in order 


to prevent a further recurrence of the in-flight failure. 


II. GENERAL AERODYNAMICS 


A. MISHAP AIRCRAFT SPECIFICATIONS 

The configurations and operational parameters of the 
mishap aircraft were similar in that both were at high gross 
weight (RAAF at 127,000 pounds, USN at 135,000 pounds) and 
operating at high airspeeds (RAAF approximately 380 knots, USN 
approximately 350 knots). Both aircraft executed a pullup 
maneuver from an altitude of less than 400 feet at the time of 
leading edge failure. The RAAF maneuver was a straight pullup 
(wings level), while the USN maneuver was a starboard rolling 
pullup (meaning that the aircraft was banked to the right as 
the pullup maneuver was executed). Interviews with some USN 
crewmembers indicate that the failure of their leading-edge 
section may have begun a few seconds earlier as the aircraft 
rolled from a right bank to wings level while inbound for the 
above stated pullup. This, they said, was reported to them by 
ground observers. Table 1 presents the relevant P-3C 
dimensional data used during this analysis while Table 2 
delineates the aircraft performance parameters. Figure 2 
shows the operating flight envelope for the aircraft in a 


flaps-up, gear-up configuration. 


TABLE 1. DIMENSIONAL DATA [Refs. 1,2,3] 


Wing 
Area, S (ft?) 1300 
Span, b (theoretical, ft) 99 
MAC, cbar (ft) 14.1 
Aspect Ratio, A Teo 
Taper Ratio, A 0.4 
Dihedral, (.25c,,, degrees) S500. 
Incidence, Root (degrees) 3.20 
Tip OLS 


Airfoil Section, Root NACA0014-1.10 40/1.051 cl.=.3,a=.8 
Tip NACA0014-1.10 40/1.051 cl.=.4,a=.8 


BE 


Straight Element OLS c 
Chord, Root (ft) 18.9 
Tip 6 
Aileron 
Area, Sa (ft?) 45.5 
Hange Line (c,,) Vet co 
Deflection Limit, Up (degrees) Zoe S 
Down FA Ole 2 
Horizontal Tail 
Area, S, (ft*) 207148 
Tae! Length (.25ecbar,, to .25cbar,, ft) 49.8 


TABLE 2. PERFORMANCE PARAMETERS [Refs. 1,2,3] 


Flight Design Gross Weight (pounds) 135000 
Load Factor (through 135,000 pounds, G) al. Omeo +35/0 
Maximum Operating Speed, Sea Level (knots) 405 
Center of Gravity Limits (135,000 lbs, %MAC) 21 Sato. 3 igx6 
Maximum Shaft Horsepower (per engine) 4600 
Maximum Lift Coefficient, C,,., (power off) 1230 
Lift Curve Slope, C,,, Tail Oft (power off) 4.84 
Tail On (power off) ~ 5 5D 
Moment Slope, C,,.,, Tail Off, c.g @ 0.25cbar, 0.19 


FLAPS UP - GEAR UP 
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ener 
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Figure 2. Operating flight envelope, clean configuration 
[Ref. 2]. 


B. PROBLEM STATEMENT 

Accurate analysis of the problem within limited time and 
budget constraints available for completion of the work served 
to restrict the options for solution methods. While 
instrumented flight and maeraloly este would probably be 
more precise, these choices were not feasible. It was 
necessary to develop an analytic method which could provide 
credible results within the academic environment. 

Determination of the loads applied and their effect upon 
the leading edge structure was accomplished by a three-step 


process. The first step was to conduct a static aeroelastic 


span-load analysis of the wing using a computer program to 
determine the section lift coefficients and structural twist 
on the P-3 wing. Next, a two-dimensional panel method was 
employed to find the pressure distribution around the leading 
edge. Finally, the forces derived from the pressure 
distribution were used as the loads applied in a finite 
element analysis computer application. mi “adareren, 
consideration was given to the aeroelastically-derived 
Spanwise wing twists which were introduced as torsion-like 


deflections in the finite element analysis. 


C. ASSUMPTIONS 

Certain assumptions were made in conducting the analysis. 
While static aeroelasticity effects upon the wing were 
accounted for, the fuselage and tail were assumed as rigid 
structures. In addition, the effect of fuselage interference 
on wing lift distribution was neglected. That is, the wing 
was taken to be fully effective in producing lift even in the 
central region where it is influenced by the fuselage. The 
unswept, straight-tapered wing with an aspect ratio of 7.5 was 
modeled structurally in the span-load analysis by an elastic 
axis. Chordwise bending of the wing was neglected. Inviscid 
solution methods were employed in the span-load and airfoil 
analysis programs. It was felt that these assumptions would 
model the flowfield without inducing an unacceptable level of 


error in the overall outcome, and that this analysis was a 


preferable method of solution to performing a computational- 
fluid-dynamics analysis where static aeroelastic influence 
could not be included. Static aeroelastic effects upon wing 
loadings were suspected of playing a large role in the problem 
because of the location of the wing’s elastic axis at a 
constant 40 percent of chord, according to available 
information. This is well aft of the 25 percent of chord 
location at which a majority of the aerodynamic loadings are 
presumed to act, creating the potential for a significant 
coupling influence by structural deflection during the 


development of lift. 





III. WING SPAN-LOAD ANALYSIS 


A. THEORY 

Because the P-3 wing (as well as nearly all others) is 
flexible, the aerodynamic-—load dds Ttesiom varies from that 
which would be seen when considering the wing as a rigid 
structure. Allowing for structural deflection of the wing 
increases the accuracy of the result in the attempt to model 
the true physics of the problem. The computer program used 
for this analysis, written in Microsoft Quick BASIC®, was 
based upon the work of Schmidt [Ref. 4]; c.£. Appendix A. The 
‘basic concepts employed in the program are presented below. 

Defining a set of linear, simultaneous equations for the 
span-load solution on a wing at a specified angle of attack 
starts with the following: 

€431(2/q)1 + Aya(2/q)2q +--- + ARn(2L/7O)y = ay 


where the wing is cut into a series of spanwise stations and 


@ a;. = aerodynamic influence coefficient of (2/q) at 


2) ; ’ 
station "j" upon induced angle at control station 
i 
e £2 = running span-load (1lb-in™?) 
®* q = freestream dynamic pressure (lb-in™?) 
Sea, = geometric angle of attack at control station "i". 


The relationship states that the span-load-induced downwash 
velocities satisfy flow tangency at a control point. 
Expressed as a matrix equation this system becomes: 


[A]{2£/q} = {a} 


where, 
e (A}] = square (n xX n) matrix of aerodynamic influence 
coefficients (length ~) 
® {2£/q} = column (n x 1) matrix of span-load values 


(length) 
® {a} = column (n x 1) matrix of angle-of-attack input 
(radian) 

In Similar fashion, a matrix equation relating the 
structural twist at a wing station "i" to the span-loading may 
be developed from: 

Si121 + Siok, tee + Sinkn = Atte 
Or, q{s;,(2/q), + Syo(2/q)o +.2. + S3,(2/q),} = Bag; 
which may be stated to apply to all the wing stations as 
before to yield the matrix equation, 
qQ{S]{2£/q} = {dag} 


where, 


e (S] = square (n x n) matrix of structural influence 
coefficients (length-force~) 


e {4a} = column (n x 1) matrix of effective angle-of-attack 
changes due to structural twist (radian) 
Next, the two matrix equations may be combined to generate 


a Single equation for an elastic wing as follows: 


10 





(A}{2/q}p = {hp +{Aa}y 
[A]{2/d}— = {%}p + Q[S]{2/q}_ 
[tal - a(si] (2/atg = {ap 
yielding the span-load solution for the elastic wing, 
(2/ate = lta) - atsil~ar, 
where the subscript "E" is used to mean elastic and "R" to 
mean rigid. 

The concept of mathematical symmetry may be employed to 
develop equations for a symmetric and anti-symmetric load case 
by considering the wing to be divided into left and right hand 
wing panels at the spanwise centerline. This process allows 
a superposition of linear aerodynamic solutions such as the 
combination of a symmetric pullup and = anti-symmetric roll 
rate to yield the total solution for a rolling pullup. It 
also reduces computing time since the size of the matrices is 
half the original required for the total wing. In employing 
this technique it is necessary to distinguish between the 
symmetric and anti-symmetric forms of the aerodynamic 
influence coefficients. 

The approach employed for developing the aerodynamic 
influence coefficients involves use of a technique known as 
the "Modified Weissinger" approach, wherein the wing is 
divided into a series of spanwise stations with swept bound 
vortices attached at the local quarter-chord point, giving 
rise to horseshoe vortices extending downstream to infinity in 


accordance with Helmholtz’ laws. The vortex strengths are 
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determined in accordance with the Biot-Savart law, such that 
the induced downwash angle at all local three-quarter chord 
points (from the summation of all vortices) are equivalent to 
the local geometric angle of attack. The three-quarter chord 
points are termed "control points". Enforcement of the stated 
boundary conditions replicates the occurrence of tangential 
flow over the surface of the wing. Application of the Biot- 
Savart and Kutta-Jukowski laws to the geometric relationships 
of the horseshoe vortices and control points results in 
development of the symmetric and anti-symmetric aerodynamic 
influence coefficient matrices, 

[Ag] =| (Alen + CAlue 

[A,] L (Alzs 7 A na 


where the subscripts "RH" and "LH" denote the right and left 


hand wing panels, respectively. These influence coefficient 
matrices may next be substituted in the previously-developed 


equations to yield, 


[A,]{£/q} 
[A,]{2£/q} = {ag} 


Subsonic compressibility effects were included in the 


1A, } 


development of the aerodynamic influence coefficients using 
the Prandtl-Glauert planform distortion approach [Ref. 5]. In 
this method, the chordwise dimension of the planform is 


increased according to the relationship, 
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The "Modified Weissenger" approach was altered to a panel 
form for this analysis in that the wing half was divided into 
five chordwise and ten spanwise stations, requiring 
manipulation of (50 x 
5a) matrices. A 


representative sketch of WING DIVIDED 


INTO PANELS 
the method employed is 


shown in Figure 3. The 
vector U in the figure 
represents the free 


Stream velocity while [I 
PANEL (I,J) 


: i A BOUNDARY 

1s the circulation 

strength of the vortex MYT HORSE SHOE 
VORTEY (1,3) 

filament. 


O8 PANEL CONTROL 


Development of the 


structural influence 





a Figure 3. Wing panel model [Ref 4]. 
coefficient matrix in 


the program is accomplished through determination of the 
moments exerted about the longitudinal (X) and lateral (Y) 
axes of the wing at points along the wing elastic axis as 
shown in Figure 4. These points correspond to the mid-span 
locations of the individual panels, along the elastic axis. 


These moments are found by solving the matrix equations: 
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{M, }=q(B/N)[t]{L/q} 
{M, }=q(B/N)*(m] {2/q} 

MOMENTS My. € My; 
Where B is the span, N ar ee 


STATION J 


is the number ops 


0O.cS Cvuoeod 


spanwise local stations, 
[t] is a torsional 


matrix and [mj is a 


ELASTIC AXIS 





bending moment matrix. 
Figure 4. Wing moments [Ref. 6] 
From these equations, a 
coordinate axis rotation can be performed as follows: 
{T} = cos, {My} + sinA,{M,} 
(wy =- sind, {My} + cosA,{M,} 
Next, torsional and bending stiffness (GJ and EI) values along 
the elastic axis are employed to determine the angular 


deflections resulting from the torsion and bending moments 


from the equations, 


Y Y 
cosa, cosa, 
O.(y)= | aes ©, (m) = | aoa” 


These are then discretized as, 


Bo MS), MZ) 
2NcosA, J2T+1 EI(J) EI(l) 





AO, (I ) a 


N 
B iy ci oma (By. 


AB (2) * ere ee GJI(J) GI(I)? 
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or expressed in matrix form, 


Be .:B il _ B 1 
{0,}= (ull a HT 1o,,) (u]( = Jim 


2NcosA, 2NcosA, 


The angle-of-attack column vector may then be found as, 


{Aa} = cosA,{0,} - sinA,{9,} 
which is also: {Aa,} = q{(S]{1/q}. 


The preceding equations may then be tied together to form the 
equation for the structural influence matrix, 

ic = (B/2N)2[uj[1/Ga]l cosA[t] + sinA.(B/2N)[m]] 

- tanA,(B/2N)2[uJ[1/EI]l sinA[t] = cosA, (B/2N) [m]]| 


B. USAGE OF THE SPAN-LOAD PROGRAM 
Application of the wing span-load program required the 
determination of the following influences: 
e Additional loading distribution due to wing angle of 
attack 
e Built-in geometric twist 
e Airfoil camber distribution 


e Dead-weight induced wing twists due to propulsion system 
weight 


e Aileron float angle 
e Propeller slipstream effects 
e Aileron control deflection 


e Roll helix angle 


Ike 


The last two influences involved anti-symmetric wing span-load 
solutions. The total wing span-load distribution, which 
provided a measure of wing lift coefficient (C,), moment 
coefficient (C,4,.) and section lift coefficient (C,), was 
obtained by an appropriate linear combination of the above 
influences. These influences were incorporated in the program 
as discrete angle of attack adjustments at the control points. 
Finally, a specified airplane lift coefficient required an 
estimate of the tail 1ift contribution before the 
representative wing C, could be estimated. The tall Iie 
contribution was based upon trimming the P-3 airplane for a 
specified flight condition and the assumption that the body 
and horizontal tail behaved approximately as rigid structures. 
1. Additional Loading 

Static aeroelastic effects upon wing span-loading due 
to geometric angle of attack were incorporated in the program 
as a selectable input from the operator at the computer 
terminal. Calculation of the appropriate angle of attack for 
a given flight condition was based on the fundamental equation 

Cr = Cro + Cred, 

with appropriate modification for tail lift contribution as 
delineated in subsection nine of this Chapter. 

Early in the analysis, it was found that static 
aeroelastic effects had a dramatic impact on additional 


loading. Solutions were obtained using the span-load program 
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for a wide range of dynamic pressures. The q equals zero 
solution corresponded to a rigid-wing case. A dynamic 
pressure of 3.7 psi corresponded too the aircraft operating at 
a Mach number (M) of 0.6 at sea level. The variation of wing- 
alone C,, with dynamic pressure, a esgs in Figure 5, indicates 
a 41 percent increase due to static aeroelastic influences at 
the sea level flight condition. The increase in 1lift-curve 
slope is associated with a spanwise variation of structural 
twist as shown in Figure 6. At a Mach number of 0.6 for sea- 
level flight, each degree of geometric angle-of-attack input 
at the wing root results in 1.95 degrees of a at the tip, with 
the added 0.95 degrees being due to the structural twist 
component in the a direction. The static aeroelastic 


influence upon the wing aerodynamic center was determined as 


| Sea Level | 
S 


per rad. 


Ae ola = 401 
0.25c Sweep = —1.3 deg 
e.g. af 0.4 ¢ 
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Figure 5. Aeroelastic variation of lift-curve slope 
(Ref. 6] 


iy 


being negligible, a result which may be attributed to the wing 


0.25 chord line having a small sweep angle. 
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Figure 6. Aeroelastic variation of twist 
[Ref. 6] 


2. Built-in Geometric Twist 

The effect of 2.5 degrees of washout was included in 
the program by linearly varying the local (panel) angle of 
attack moving outward from zero at the root to -2.5 degrees at 
the tip. Figure 7 depicts the effect of washout on section 
lift and twist distribution for the rigid and elastic P-3 
wing cases. The magnitude of the negative section C, values 
is seen, in Figure 7, to increase by 21 percent while the wing 


tip experiences a 1.5-degree negative twist due to aeroelastic 


influences. 
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Figure 7. Effect of built-in geometric twist on section 
lift coefficient and twist distribution. 


3. Engine Dead-weight Twist 

A separate version of the span-load program was 
developed and run to determine the wing twist induced per G 
due to the effect of engine dead-weight moment. Output 
information was generated in terms of discretized angle-of- 
attack adjustments and included in the basic’ span-load 
program. An engine and propeller assembly weight of 3974 
pounds was assumed to act at the (x,y) coordinates (coordinate 
system as shown in Figures 3 and 4) of (-51", 187") and (-45", 
357") for the inboard and outboard engines, respectively. 
Figure 8 shows the effect of engine dead-weight twist on the 


elastic twist distribution at three G’s and 405 knots. The 


1g 


wing is twisted to a maximum of -4.5 degrees at the tip under 


this load condition. 


N) 
Uy 


—+— Rigid Twist 


a 
i 
L 
fom 
2 
a 
=. 
2 
3 
ho 


i 
ae) 


—O— Elastic Twist 


Wing Station (2y/b) 





Figure 8. Effect of engine dead-weight on structural 
twist. 


4. Camber 

Although listed in Table 1 as a NACA 0014 which tapers 
to a NACA 0012, the airfoil used for the P-3 wing is not a 
symmetric section, as indicated by available data [Ref. 8, p. 
2-470]. It is, in fact, a hybrid with considerable camber. 
The wing camber has an effect on aeroelastic behavior and was 
included in the analysis by determination of camber-line slope 
at the five chordwise control points which make up the wing 
model and entering the negative of this value as an adjustment 


to the station angle of attack. 
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5. Aileron Float 
The span-load analysis included the effect of aileron 
float as given in a graph of aileron angle versus airspeed in 
a Lockheed report [Ref. 9, p. 123]. An equation curve fit for 
the 50-pound control-wheel-force curve was developed and 
incorporated in the program. Deflected aileron surface area 
was matched by a selection of panels on the wing model for 
varying degrees of deflection in order to closely emulate the 
aircraft control-surface deflection and aeroelastic effect. 
6. Aileron Deflection Angle 
Using the same process as that used for the aileron 
float angle, the effect oof aileron deflection for 
consideration of anti-symmetric solutions to emulate roll 
maneuvering was also included. A curve fit to the data given 
for available 50-pound control-wheel-force deflection in the 
Same report [Ref. 9, p. 123] was used to generate the 
deflection angles. 
7. Roll Helix Angle 
Data [Ref. 9, p. 118] for available tip helix-angle 
(pb/2V) variation with airspeed were curve fitted and 
incorporated as another control point angle-of-attack 
variation in the program. A roll helix angle of 2.63 degrees 
provided roll moment equilibrium with the available 9.3 
degrees of aileron deflection at 275 knots equivalent 


airspeed, at a control wheel force of 50 pounds. At 405 
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knots, the roll helix angle was 1.03 degrees for roll 
equilibrium with an aileron deflection of 8.2 degrees. The 
merging of the anti-symmetric input due to aileron deflection 
with the symmetric pullup solution allowed a comprehensive 
analysis of the static aeroelastic effect of a rolling pullout 
maneuver as described in part A of this Chapter. An example 
of the individual and combined effects of these components on 
wing section lift coefficient (C,) distribution is shown in 
Figure 9, where the outer corner of the operating envelope for 


roll maneuvering (405 knots, 2.4 G’s) was examined. 
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Figure 9. Rolling pullup contributions to section lift 
distribution. 


8. Propeller Slipstream Effect 
The increase in dynamic pressure over the wing due to 


the propellers was calculated using momentum theory as stated 
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in Glauert [Ref. 10, p. 200] according to the following 


formula: 
T = AAV + V1)Vv1 
where, 
e T = propeller thrust (pounds) 
e A = propeller disk area (ft?) 
e v = freestream velocity (ft/sec) 


® v, = velocity increase behind the propeller, determined 
from known thrust. 
The velocity increase, v,, was converted to a dynamic pressure 
boost and incorporated in the span-load program in the area 
behind the propellers. 
9. Tail Contribution 
Horizontal tail and fuselage moment effects on wing 
span-load distribution resulted in an adjustment to the wing 
angle-of-attack value used as input to the program. This 
adjustment was calculated according to the relationship, 
Cr = Craaa + Cre + Cro 


where, 


ec, = lift coefficient required for flight condition 


© Cragq = Lift coefficient due to additional loading (due to 
angle of attack) 


Seen — tall contribution to lift coefficient 


® Cro tail-off lift coefficient at zero angle of attack 


(due to camber, wing twist and dead-weight twist) 
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and, 


Cie = “MM, (char 2, ) 


where, 


e AC. = tail-moment coefficient 
e cbar = mean aerodynamic chord of the wing 


ef, = distance from .25cbar,, to »25cbar, 


Additionally, at airframe trim (C,cg, = 0), 


Crt Cmo + Cmercnada 


where, 


© C,9 = moment coefficient at zero angle of attack 


® Cie, = tall-off airplane dcC,/dC,, c.g. = 0.25cbar 


The above was assembled and solved for C,.4, to give, 


cbar 
Cos, Sigal T, Cmo 





Cradd= 
4+ Char 


CG 
7 - mCL 





The value of tail-off C,,, for the aircraft was found from 
wind tunnel data to be approximately 0.19 [Ref. 3, p. 20]. 
The value of C,.4, for any given flight condition was then used 
to solve for the angle of attack in the span-load program 


according to the formula: 


a = Craaa/ra 
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The value of C,, was available as an output, for the elastic 
wing, from the program for any given airspeed by entering one 
radian for the angle of attack input. 

Becdisacct longer 0503 CoC jadue tomfuselage effect 
was made after finding that the value C,, in the wind tunnel 
data and that from the program differed by this amount. This 
correction was verified by calculations made in accordance 
with Etkin [Ref. 7, p. 334] concerning the effect of body and 


engine nacelles on neutral point location. 


C. EFFECTS NOT INCLUDED 

Two other factors were considered as possible contributors 
to the static aeroelastic problem, but were found to have no 
Significant impact and therefore not included in the span-load 
program. These factors were propeller gyroscopic effects and 
the effect of moments due to wing fuel. 

1. Propeller Gyroscopic Moment 

This phenomenon was investigated using the following 


formulation from [Ref. 11]: 


e M = moment due to gyroscopic precession (lb-ft) 
p = polar moment of inertia of each blade (lb-ft sec“) 


® ® = rotation rate of the propeller (rad-sec’+) 
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® 2 = pitch rate of the aircraft (rad-sec™1) 
e m = mass of propeller blade (slugs) 
e k = radius of gyration (ft) 


The analysis showed that, allowing for a five-degree-per- 
second aircraft pitch rate, the moment developed by each 
propeller was 1225 lb-ft in a counterclockwise direction as 
viewed from above the propeller. This moment was considered 
insignificant since it is applied in the lateral plane and 
does not influence the static aeroelastic span-load solution. 
2. Wing Fuel Moment 

Analysis of the fuel tank geometry and location 
revealed that the center of mass of the fuel (with full wing 
tanks) lies nearly coincident with the elastic axis. Any 


torsional moment derived from this source would be negligible. 


D. APPLICATION 

Linear summation of the contributing factors addressed 
earlier in this Chapter resulted in a prediction of the Loan 
span-load distribution for the P-3 wing under any given flight 
condition. The primary flight conditions of concern in this 
analysis were those encountered during the aircraft mishaps as 
discussed in the introduction. The basic premise applied was 
to examine limit loads at the edge of the operating envelope 
and then expand the analysis to regions outside the envelope, 


at the ultimate load condition. In addition it was decided to 
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first examine the loads encountered in a symmetric pullup of 
three G’s at 275 knots. This speed was chosen because of an 
observed higher load on the leading edge than determined at 
the upper right corner of the flight envelope, at 405 knots. 
(A further discussion of this loading phenomenon will be 
presented in Chapter IV.) Following that, the target flight 
condition was extended to 4.5 G’s at 325 knots, which 
represented an approximate extension of the envelope uSing a 


value for C of 1.3. Next, the effect of anti-symmetric 


Lmax 
loading in the form of a starboard and then a port rolling 
pullup were considered in order to assess the contributions of 
aileron deflection and roll helix angle. Rolling pullup 
analyses were done at 2.4 G’s and 275 knots to remain inside 
the envelope and maintain some congruity with the symmetric 
loading case. Presented in this section are plots of the 
impact of some of the various flight maneuvers on section lift 
coefficient and twist distribution, using the fully-developed, 
tail-on solution with all contributing factors included. 
Figures 10 and 11 show the 275-knot, 3-G symmetric condition, 
and compare section lift coefficient and structural twist for 
the rigid-wing and elastic-wing cases. In Figures 12 and 13, 
the same distributions are depicted for the 275-knot, 2.4-G 


rolling pullup load condition. 
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Figure 10. Lift coefficient distribution at 275 knots, 
3-G, symmetric pullup. 


——*— Rigid 


—— Elastic 


Twist (deg) 


efeeeeeq@ee@ee@ @Fwewaswe @M@ewsz_ BSieeeees =a = 


0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 
Wing Station (2y/b) 


Figure 11. Twist distribution at 275 knots, 3-G, 
symmetric pullup. 
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IV. TWO-DIMENSIONAL PRESSURE ANALYSIS 


A. THEORY 

Like the static aeroelastic span-load analysis, the two- 
dimensional pressure analysis employed in this thesis involves 
the application of linear superposition. This is a 
consequence of the pressure analysis being based upon a 
solution to the LaPlace equation, a linear, homogeneous 
second-order partial differential equation. Linearity allows 
the problem to be subdivided into three separate elements 
(which will be described later in this Chapter) and added 
together. The actual application was based upon a panel 
method similar to the technique used in the preceding Chapter 
except that now, instead of dividing a wing planform into 
chordwise panels, a two dimensional (x,y) airfoil is divided 
into panels along the perimeter of its surface. Camber is 
defined in the airfoil shape, instead of being added on as a 
discretized variation of angle of attack, as before. The 
analysis now concerns a vertical plane or cross-section. 
Steady (no variation in the flow field with time), inviscid, 
incompressible flow is assumed to exist. The panel method and 
formulation employed is documented in a Naval Postgraduate 


School thesis by Teng [Ref. 12]. 


30 


1. Coordinate Axis and Panel Numbering System 

The airfoil is considered to be fixed in an (x,y) 
coordinate system with its origin at the point of intersection 
of the chord line with the leading edge. The positive x axis 
points aft to the trailing edge while positive y is up. The 
panels which make up the airfoil surface are of varying 
lengths, depending mostly on the radius of curvature, and are 
numbered from 1 through n starting at the lower surface of the 
trailing edge and proceeding clockwise to the upper surface at 
the trailing edge. Delineating the end points of these panels 
are nodes which begin with the number 1 at the trailing edge 
and proceed Alene the same numbering path as the panels. The 
trailing edge point is counted twice, giving n+l nodes in all. 

2. Flow Formulation 

Consider some panel j on the surface of the airfoil. 
On this surface there exists a pair of singularity 
distributions, known as a source distribution q; and a 
Memeteity distribution jy. The strength of the source 
distribution varies from panel to panel while the vortex 
Strength is the same for all panels. These singularity 
distributions satisfy LaPlace’s equation and the far field 
boundary condition. 

Applying superposition, the overall flow field is 
considered to be made up of three individual flows and is 


represented by the equation, 


5) JI 


@= $+ 5 + by 
where ¢, is the potential of the freestream flow, 
do = V(X cosa + y Sina) 
¢?, is the velocity potential of the source distribution of 


strength q(s) per unit length (s) and is calculated by, 
s 
$= [| US) 1n(r)ds 


where (r) is the radial distance from some point at which a 
source and vortex flow exist, to the midpoint of the panel in 
consideration. In addition, ¢, is the velocity potential of 
a vorticity distribution of strength X(s) per unit length and 


is given One 
> 21 


where 0 is the angle formed by a line drawn along the radial 
distance (r) and the panel in question. 

Each of the preceding equations is integrated along 
the straight line which makes up each of the panels, where q; 
and Aare constant. The individual effects are then summed to 
give the total effect of the sources and vortices from all 


panels, as given in the equation, 


Nn 
; 5 a4 
S=V ae -_t 
a(xcosatysina) + )/ | , [ sz 1n(r) 5-§]ds 
J=1 panél(j) 
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The calculation of ® requires solution of the (n+1) 
unknowns, q; (ja 1,2) anche. This is accomplished 
numerically in the computer program. Once %@ is known, the 
velocity can be found by taking the gradient (V) of @. The 
total velocity vector is found as, 

Veotal = V@ = Von + V(os + $y) 

Next, the coefficient of pressure is found from the 

Bernoulli equation in the incompressible form, 
al ia anes )2 
3. Boundary Conditions 

Both the condition of flow tangency at the surface and 
the Kutta trailing edge condition must be satisfied as 
boundary conditions. As in the span-load analysis, control 
points are designated at which flow tangency is satisfied, 
except here the control point is taken as the mid point of the 
panel. It is stipulated that each control point will have 
tangential velocity, > ., but that all normal velocities, 
(V");, will be exactly zero. 

The Kutta condition requires that the pressures on the 
upper and lower surface at the trailing edge be equal. Using 
Bernoulli’s equation for steady potential flow, this state of 
pressure equilibrium is found to exist when the tangential 


velocities in the downstream direction are equal at the upper 
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and lower trailing-edge panels. In equation form this is 
written, 
Cy Sy 

The task then becomes one of using the boundary conditions to 
solve for the (n+1) unknowns. 

4. Influence Coefficients 

The concept of influence coefficients is again 
employed in this portion of the analysis, as it was in the 
span-load analysis. Here, the influence coefficients take 
the form of induced normal and tangential velocities at the 
control point of a given panel. These velocities are induced 
by the source and vorticity distributions of the other panels, 
and it is from this influence that they receive their 
designations: 
e A"; 5 = normal velocity induced at the i* * panel control 
point by the source distribution on the j* panel. 


e A‘,, = tangential velocity induced at the it? panel control 
point by the source distribution on the j*" panel. 


e Bee normal velocity induced at the i&® panel control 
point by the vorticity distribution on the j*? panel. 


e Bt, , = tangential velocity induced at the i*? panel control 
point by the vorticity distribution on the j* n panel 
These influence coefficients are calculated through 
application of the geometric relationships which exist between 
the panels in conjunction with the formulas for the source and 


vorticity velocity potentials as given above. 
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5. Numerical Solution Method 
Using the influence coefficients, the boundary 
conditions can now be employed to write (nt+1) equations which 
may be solved in matrix form for the (n+1) unknowns. The set 
of n equations comes from enforcement of the flow tangency 


boundary condition in the form, 


Nn 


Nn 
dD [A* 5595] +y B® ,3+V,sin(a-8;) =0 
j=1 j=1 | 


Next may be written the enforcement of the Kutta boundary 


condition as, 


Nn fl 7} Nn 
eee) 1590,)-Y ), BX1;-Vocos(a-0,)=) [A*,7g5] ty}, Bn7t¥_cos(a-6)) 
j=l jal jal Jae 


The negative signs on the left side of the equation are due to 
the defined orientation of the tangential velocities as 
positive in the downstream direction. 

These equations may then be expressed in matrix form 
with the (n+1) unknowns (il.e., q; (Gat 2,6.) and Y) 
arranged as a column (nxl1) matrix multiplied by the Aji4 psi 
(nt+1 x nt+1) influence coefficient matrix and set equal to a 
Bu+, column matrix. From this point, a Gaussian Elimination 
numerical technique may be employed to solve for the (n+1) 


unknowns [{Ref. 13]. 
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6. Velocity and Pressure Distribution 
Once the q, and y are found, the tangential velocities 


at the control points can be solved for according to the 


equation: 
N nN 
Vectal=», [A*ssq5]+*¥ > B®, ,+V,cos(a-6;) 
oul oa 
wher@G wk = 17275... n- From this, the individual pressure 


coefficients may be solved using the equatioii, 

(Co): Sela 2) | 
At this point, the forces at work on the airfoil may be found 
by first integrating forces in the airfoil coordinate system 


as follows, 


Nl 
Cy de (Cp) (541 -X3) 


Nn 
C= yo (CriGg taieD 
1=1 


Performing a coordinate axis rotation to re-align with that of 
the freestream yields the lift coefficient, | 

Cp, = Cy cosa - C, Sina 
For this application, the values of C, an C, are found in 
discretized form at each panel as (C,); and (Cy)ay i = 
1,2,..-.,n, and then multiplied by the chord length and 
freestream dynamic pressure to arrive at the normal and 


chordwise force exerted at each panel. From this, the 
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leading-edge panels are selected and their forces collected 
for application as point loads in the finite element analysis 


portion. 


B. APPLICATION 
1. P-3 Airfoil Section 

As mentioned in Chapter III, the airfoil shape used in 
this application is as delineated by ([Ref. 8]. Surface 
coordinate locations were solved using the tables and 
equations provided therein relative to a "wing reference 
plane" which appeared to correspond to a water line (i.e. the 
angle of incidence at the root was included in the 
definition). These coordinates were then rotated to an (x,y) 
coordinate system aligned with the chord of the airfoil as 
required in this panel method. The result was an airfoil 
consisting of some 48 panels, which was increased to 52 panels 
(53 node points) due to observed roughness of the leading edge 
shape when plotted. This smoothing of the leading edge shape 
was achieved by applying the specified leading edge radius to 
create intermediate node points. The basic shape of the 
airfoil, though not precisely to scale, is depicted in Figure 
14. The locations of the node points are also shown. 

2. Program Inputs and Outputs 

The desired output from the two-dimensional panel 

method program was a collection of forces and boundary 


constraints to be applied at node locations in the finite 


on 





Figure 14. P-3 airfoil showing panel nodes. (Not to 
scale) 


element analysis phase. In order to achieve this, the program 
was altered considerably from the initial state as outlined 
earlier in this Chapter. The final FORTRAN code is available 
in Appendix B. Since it was necessary to develop loads to be 
applied to a three~dimensional model, the program was set up 
to iteratively compute the force distribution at a series of 
two-dimensional airfoil sections which ranged in location from 
the inboard end of the finite element model at wing station 
(WS) 256 to the outboard end at WS 320. (These locations 
correspond to the outboard 64 inches of the wing center 
section leading edge, located between the nacelles.) This 
section force distribution was scaled up from a chord- 
normalized airfoil shape to a full-sized airfoil as described 


above, then multiplied by a scaling factor equal to the 
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distance between the section locations. Section locations 
consisted of 24 wing-station positions along the leading-edge 
model, including the nine ribs and 15 intermediate points as 
determined by finite element model node locations. 

The first of the inputs to the program consisted of 
the airfoil geometry definition as described above. Next, a 
tabular file of wing station locations was read in together 
with.the section lift coefficient for that spanwise location 
as found by the span-load program. These section lift 
Ceeaticients, initiaally found for the .35, .45 and .55n 
locations (where n= 2y/b) were curve-fitted using the Cricket 
Graph” plotting software in order to achieve a high degree of 
accuracy in determining the individual C)’s at stations which 
were no more than three inches apart. Also included in this 
input file was a list of spanwise multiplication factors for 
scaling up the load as described above. Additional input 
files consisted of the finite element node numbers which were 
matched with their respective x and y direction loads in the 
program. Read in from the terminal were the airspeed under 
consideration and the twist angles of the wing box as 
determined from the span-load program. 

Output consisted primarily of the load file which included 
not only the loads at the finite element node points, but also 
the constraints and twist displacements for the upper flange 
and lower hinge node points of the leading-edge finite element 


model. The twist displacements were calculated within the 
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panel program using formulas based on geometric considerations 
of the height of the wing spar at the inboard and outboard 
ends of the leading-edge segment, and the net twist 
displacement of the front spar from the inboard to the 
outboard end of the model. That is, the finite element model 
was assumed to have undergone a rigid-body rotation to the 
degree of twist which was found to exist at the inboard end 
(WS 256). Therefore, twist displacements at the inboard end 
were set to zero, followed by application of the subsequent 
net twist distribution that occurred at the other wing 
stations while proceeding outboard to WS 320. This net twist 
distribution was equal to the difference between the outboard 
and inboard twist amounts, applied linearly over the 24 
Stations. This twist amounted to approximately 0.3 degrees in 
the 275-knot, 3-G symmetric pullup. The displacements 
generated for application at the node points were in the 
longitudinal (x) and vertical (z) directions on the three- 
dimensional model. Twist rotation of the front spar was taken 
to be about its vertical mid point, and the structural wing 
box was assumed to have no chordwise distortion — rotated 
about the .40c elastic axis location. Other output took the 
form of files to examine Cy distributions and to tabulate 
loads by wing station and two-dimensional node point for 
verification of the load file. In addition, output was 
generated which approximated the total normal and chordwise 


loads applied to the entire wing leading-edge center section 
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located between the engine nacelles. This estimation was 
accomplished by computing the load on the leading-edge portion 
of the airfoil at WS 274 and multiplying by the length of the 
leading-edge segment (92 inches). 
3. Program Operation 

After reading in the coordinate information for the 
airfoil along with the wing station, C, and _ load 
multiplication factor, the FORTRAN program calculated the 
chord length at the particular station based upon the 0.40 
taper ratio. It also determined the thickness fraction of the 
airfoil section at that station by assuming a linear taper 
from 14% maximum thickness at the root to 12% at the wing tip. 
This thickness factor was then used to recalculate the y 
coordinate position of each node point to redefine the shape 
of the airfoil. An initial angle of attack of one degree was 
set and the process described in the theory section of this 


Chapter took place, wherein the C,, C and Cy, were calculated 


y! 
for that angle of attack. This C, value was then compared to 
that required (as input with the wing station), and an 
iterative cycle commenced in which the angle of attack was 
varied up or down by an amount based on the product of the Cy, 
deviation multiplied by a preset angular value. An accuracy 
test of .0001 was applied to reach an acceptable value for C,, 


at which time the process started over with the next wing 


station. Forces in the x and y direction were matched with 
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the appropriate finite element node points and written to the 
load file for each iterative cycle. After all loads were 
calculated and stored, the twist displacements and zero 
boundary constraints were calculated and appended to the load 
file. 
4. Verification of the Program 

The accuracy of the two-dimensional panel method was 
verified by comparison with published empirical data for 
tangential velocity and/or pressure coefficient distributions 
for the NACA 0012 and Eppler E64 airfoils before its use in 
this application. Results were nearly identical to the 
published data with only a small deviation seen near the 
trailing edge of the program tangential velocity distribution 
for the Eppler airfoil. No difference from the NACA 0012 Cy 
data could be identified. 

5. Flight Regime Selection 

It was found in the course of running the program at 
various airspeeds and angles of attack that the highest loads 
on the leading-edge segment were generated at slower airspeeds 
and higher angles of attack as a constant G load was 
maintained on the aircraft. This result seemed contrary to 
conventional opinion that the highest loads would most likely 
occur at or near the high speed end of the operating envelope. 
A brief study was undertaken to determine the cause of this 


Phenomenon and a hypothesis is given here. 
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a. Method Employed 

In an effort to study the effect of dynamic 
pressure (airspeed) and angle of attack on the leading edge of 
an airfoil, the variation of net Cy, distribution on the 
leading edge at various angles of attack was first examined. 
These C, values were the numerical sum of the difference 
between the upper and lower Cy’S, integrated over the 
chordwise distances occupied by their respective panels. 
These data obtained were then curve fitted with a third order 
polynomial and used in a spreadsheet to calculate the Cy 
distribution at varying angles of attack. These angles of 
attack were generated by varying the airspeed from 275 to 425 
knots, calculating dynamic pressure (q) for a 3-G wing loading 
from the Bernoulli equation, and then converting these q’s to 
angles of attack required using the basic lift formula, 
altered by the equation, 

Cr = Cro + Cree 
to give, 
a = 3W/C,.q58 and a = 3W/C,.dS8 - Cro/Cre 

for the symmetrical NACA 0012 and cambered P-3 airfoil, 
respectively. These angles of attack were then used as inputs 
to the polynomial curve fits, from which a corresponding set 
of C, values were calculated. Next the product of Cy and q 
were found, to give the pressure acting on the leading edge, 


corresponding to a matched set of q and angle of attack. This 
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information was then plotted as seen below in Figures 15 and 
16. 
b. Symmetric Airfoil 
Figure 15 shows the distribution of the various 
parameters on a relative scale as airspeed increases for the 
case of a symmetric airfoil. Note that the pressure acting on 
the leading edge is nearly constant, showing only a slight 


increase with increasing airspeed. 
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Figure 15. Pressure analysis of 0012 leading edge. 


c. Cambered P-3 Airfoil 
Next the same information is plotted for the P-3 
airfoil in Figure 16. Note that the pressure on the leading 
edge undergoes a marked decline as airspeed increases (AOA 


decreases). 
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Figure 16. Pressure analysis of P-3 leading edge. 


d. Analysis 

Because of the apparent drop in pressure on the 
cambered airfoil, it was concluded that the effect of camber 
was to shift the loading of the wing in a way that caused 
angle of attack to become the dominant influence rather than 
dynamic pressure. It was also observed that if the lift 
equation applied in the program was that of 2 cambered airfoil 
(i.e. Cyy had some positive value as the lift curve was 
displaced upward), this loading phenomenon was present. If 
Cog were zero, the load on the leading edge was approximately 
the same at various angles of attack and airspeeds. Because 
of this observation, it was decided to select the 275-knot, 3- 
G flight position as the maximum loading position in the 


normal operating envelope. 
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V. FINITE ELEMENT ANALYSIS 


A. THEORY 

The method of finite element analysis is a means of 
Simulating the structural behavior of a continuous physical 
system by a discretized representation of that system. 
Structures are represented by discrete node points which are 
connected by structural elements. The nodes form a grid which 
details the general shape of the structure while the elements, 
although they appear as only lines, are mathematically given 
the physical properties of the portion of the structure which 
they are there to represent. A physical structure is thus 
transformed into a mathematical representation for the purpose 
of analyzing some behavior of the structure. This analysis 
may be in the area of dynamic response, heat transfer, or, as 
is the case here, static loading response. This method of 
analysis is widely proven to be highly accurate and has been 
used in many engineering fields, including aerospace, 
automotive, civil and mechanical applications. With the 
increased capability, speed and data storage capacity of 
microcomputers, this analysis technique feene longer limited 


to mainframe applications, as was the case a few years ago. 
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1. Accuracy 

The accepted rule of thumb in finite element modeling 
is that the use of more node points results in a more accurate 
solution. Convergence tables have been developed which show 
that the use of fewer nodes increases the stiffness of the 
model . The main drawback to using a large number of nodes is 
that it greatly increases computation time and requires larger 
amounts of storage space than a model of the same structure 
using fewer nodes [Ref. 14]. The finite element analysis 
software used for this thesis is called MSC/pal 2° and is a 
product of the MacNeal-Schwendler Corporation, the company 
which also creates the highly respected NASTRAN® finite 
element application for VAX/VMS work stations and mainframe 
operations. The accuracy of MSC/pal 2° has been tested and 
documented by the manufacturer [Ref. 15]. In addition, the 
manufacturer recommends simple hand calculations to verify 
that results obtained for a given model are reasonable (i.e. 
within the same order of magnitude). This was done by 
calculating simple beam bending stress with constant area 
cross sections approximating the rib legs of the model. 
Results established that the finite element model produced a 


solution which was well within expected norms. 
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2. Equations 

The number of equations to be solved in the finite 
element analysis is equal to the number of degrees of freedom 
in the model. Each node point has six degrees of freedom: 
three translational (in each of the (x,y,z) coordinate 
directions, and three rotational (about each of the coordinate 
axes). Stiffness equations are generated for the stiffness of 
each connecting element, based on the specified material 
properties (Young’s modulus, shear modulus, mass density, 
tensile yield stress are specified) and the geometric 
configuration of the element. Elements may take the form of 
beams, triangles, quadrilaterals and others. These nodal 
stiffnesses are combined to form a system stiffness matrix, 
(K], of size (NxN) where N is the number of equations. 
Degrees of freedom may be eliminated by setting them to zero 
in the model definition phase (a way of applying boundary 
constraints) or fully retained as was done in this 
application. (Here, boundary constraints were applied in the 
load as discussed in Chapter IV. This application allows 
variation of the boundary displacements from one load to the 
next and allows recovery of reaction forces at the constrained 
nodes.) Once the stiffness matrix has been formed, the static 


analysis may be performed according to the following equation: 


[K]{U} = {F} 
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where, 


e {F} column vector of applied loads (Nx1) 


e {U} resultant column vector of nodal displacements 


(Nx1) 
Gaussian elimination is employed to solve the matrix equation. 
In large models, as is the case here, matrix partitioning 
takes place prior to solution. 

Once the displacements are known, stress-strain 
relationships are employed to compute stress values throughout 
the structure. These stresses are available to the user in 
the form of major and minor principal axis stresses, Von Mises 
stress concentrations and maximum shear. stresses. ite! 
addition, output of displacements and rotations’ are 


accessible. 


B. STRUCTURAL REPRESENTATION 
1. The Leading Edge Structure 

The leading edge segment considered for this analysis 
was the port wing, center section, located between the number 
three and four engine nacelles. Figure 17 shows a cutaway 
view of the structure. The segment is composed of 12 vertical 
ribs supporting a double (inner, outer) skin. The outer skin 
is .040 inches thick and the inner is a stamped corrugation of 
.016 inches. Assembly of the outer and inner skins provides 


a series of ducts approximately .25 inches in height and two 
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Figure 17. P-3 leading edge structure [Ref 16]. 


inches wide each. These ducts run longitudinally along the 


inner surface of the leading edge and are continuous over the 
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Spanwise length of the leading edge. The purpose of this 
series of ducts is to provide a channel through which bleed 
air may travel to heat the leading edge. The tube seen 
extending through the structure delivers the bleed air to the 
ducts. The leading edge segment is secured to the front spar 
of the wing by a full-length piano hinge at the bottom edge 
and a screwed-down spar cap flange at the top. This 
arrangement provides access to the area for maintenance 
functions. 
2. The Model 

The basic finite element model used in this analysis 
was obtained through translation of a NASTRAN® finite element 
model using a function in the MSC/pal 2° application called 
NASPAL®. NASPAL® reads the NASTRAN® text file for the model 
and rewrites it in the format used by MSC/pal 2”. The 
NASTRAN” model file was obtained from Aerostructures, Inc. 
through contact with NAVAIR’s AIR-530 office. Because the 
NASTRAN’ model consisted of 2749 nodes, it was necessary to 
reduce the model size to meet the MSC/pal 2° limitation of 
2000 nodes. This reduction was done by entering a set of 
geometric coordinates during the NASPAL° translation and 
instructing the translator to consider ony the outboard 64 
inches of the model. In effect, the leading edge section was 
severed between WS 247 and WS 256, or between the third and 


fourth ribs from the left end as shown in Figure 17. The 
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remaining portion of the model was translated from the fully 
defined model. This approach was considered to provide a more 
accurate solution than increasing the spacing between nodes, 
keeping the highest available 
level of detail in the model 
(again, more nodes '-= mean 
better accuracy for the same 
structure). In doing so, it 
was necessary to modify the 
inboard end rib of the 
structure (WS 256) since the 
end ribs were constructed 
differently than the 
intermediate ribs. The end 
ribs are closed in the front 
and have single, riveted 
flanges on their interior 
(Figure 18, bottom), while 
the intermediate ribs (Figure 
18, top) have an open front 
to allow access for the bleed 
air "pump cap" assembly which 
leads to the double skin 
described earlier. The 


intermediate ribs also have 





Figure 18. Intermediate and 
double flanges as_ shown. end rib detail. 
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These changes to the model were made in the rib at WS 256 (now 
the inboard end rib of the model) by adding node points and 
quadrilaterals to close the front and by rearranging the beam 
element properties within the MSC/pal® 2 model definition text 
file. In this way, the model was altered to represent a 
shortened leading edge segment with properly defined ribs. 
The model was constructed as an all-aluminum structure, and as 
stated in Chapter IV, the loads applied were generated for 
each specific wing station and scaled according to the 
Spanwise distance between the node points. In this way, the 
model detailed here received a scaled-down load for its 
scaled-down size. The upper and lower surface views of the 
model are presented in Figures 19 and 20. The end views of 
the model are shown in Figure 21. 

The lower hinge is replicated in the model by leaving 
the Y-axis rotation unrestrained. The upper flange of the 
model is secured in all six degrees of freedom. Displacements 
for front spar twist are incorporated in the upper and lower 
constraints as detailed in the previous Chapter. All of these 
boundary conditions are input through the "Displacements 
Applied" command section of the load file. 

Another difference between this model and the original 
NASTRAN code developed for NAVAIR is that the original did 
not incorporate stiffness generation in the skin of the model. 
The skin thickness of the NASTRAN® model was .056 inches, 


which is the sum of the outer and inner skin thicknesses but 
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Figure 19. Upper surface of finite element model. 
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Lower surface of finite element model. 
5 


Figure 20. 


weet 1 OU t the 
stiffness 
generation enabled 
in the model 
definition file, 
the skin would have 
no stiffness. In 
effect, it would 
ae as a non- 
loadbearing 
membrane. In the 
MSC/pal 2° model, 
stiffness 
generation was 
enabled, but when 
loaded, this 
resulted in 
excessive 
deformation of the 
skin. It was 
decided that the 


.056 inch thick 





Figure 21. End views of finite element 
model (outboard, inboard). 


skin did not accurately model the combined effect of the inner 


and outer skin combination since the corrugated inner skin 


would be much stiffer than its mere thickness (.016 inches) 


would represent. 


Although there is a variable stiffness 
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factor in the MSC/pal 2° application, available references on 
finite element analysis did not detail the calculation of this 
factor. Therefore, it was decided to increase the skin 
thickness to an equivalent thickness which would accurately 
replicate the behavior of the corrugated skin combination. 
This thickness was calculated by equating the moment of 
inertia of a box beam formed by a single corrugation and outer 
skin combination with that of a solid cross section about the 
same reference axis (at the upper surface). After allowing 
for spaces between the corrugations where the two skins are 
riveted together, an equivalent skin thickness of .1838 inches 
was determined and employed. This action reduced the 
deformation of the skin under load to reasonable norms. 

Loads were applied to the model at spanwise rows of 
selected node points which most nearly corresponded to the mid 
points of the panels in the two-dimensional panel method. 
These point loads were applied in the vertical (2) and 
horizontal (X) directions in units of pounds force. 

Material properties used to represent the 2024-T6 


aluminum structure in the construction of the model were as 


follows: 
e Young’s modulus (E) = 1.06E+07 psi [Ref. 17] 
e Shear modulus (G) = 4.0E+06 psi [Ref. 17] 
®e Poisson’s ratio (vU)= 3.25E-01 


¢ Tensile yield strength (o = 4.7E+04 psi [Ref. 18] 


ys? 


56 





where Poisson’s ratio (v) was calculated according to the 


standard formula: 


ee. ee 
2( Lave 


C. APPLICATION 

Many load conditions were examined in the course of this 
thesis. Presented here are the six load cases which give the 
best overall illustration of the observed effects of static 
aerodynamic loading, both with and without static aeroelastic 
effects included. In this way, the effect of wing box twist 
may be seen, along with the combined effect of angle of attack 
and dynamic pressure. In Table 3, the features of the six 
load cases are given. The L/R column provides a distinction 
between a wings level (L) pullup or a rolling (R) pullup. The 
SeGoet OL rolling into a turn during the application of G 
loading (as in a climbing breakaway maneuver) is not the same 
as that of rolling out of a turn during G application (as in 
rolling to wings level while pulling out of a dive). Since it 
was found that the loading effect in terms of both twist and 
(air loading was greater during the former (due to the combined 
effect of roll helix angle, aileron deflection and air load), 
the former was chosen for presentation here in the rolling 
load cases. All load cases were generated at 135,000 pounds 
gross weight except for number 4 which was done at 110,000 


pounds. In load case number 1, the effect of twist was 


3) 


eliminated from the load solution by setting the MSC/pal 2° 
displacements to zero so that the effect of wing torsional 


twist may be seen by comparison with load case 2. 


TABLE 3. LOAD CASES EMPLOYED 


Load Case Airspeed(kts) Load Factor L or R Comment 
1 275 3.0 Level No twist. 
2; 275 3.0 Level 
3 275 2.4 Roll 
4 240 Suul) Level 110,000 # 
5 350 Deg Roll 
6 325 4.5 Level 


1. Input Data 

Presented in Table 4 are the input data sets for each 
of the load cases. These are given for the three pertinent 
non-dimensional spanwise wing stations (2y/b) as generated by 
the span load program. The 24 actual wing stations (2y/b = 
0.4322 through 0.5397) employed for the two-dimensional panel 
program were solved by curve fitting the lift coefficients 
bounded by these extremes, as described in Chapter IV. The 
92-inch loads are the approximate longitudinal (positive aft) 
and vertical (positive up) total loads which would be seen by 
a complete center section leading edge segment on the 
aircraft. The total load applied to the shortened finite 


element model would be some 69.6 percent of this stated load. 
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TABLE 4. FINITE ELEMENT ANALYSIS INPUT DATA 
Load Case 2y/b (Uy Twist (deq) 92" Load (lbs) 

1 ETS hs 2 AO OECOO0 Pee = —5 157 
~45 1.2379 0.0000 FZ = 13155 
255) 1.1834 0.0000 

2 35 1 Pee PAS 0.8736 Pee — 5 1 SY 
045 ie 23 19 131221 Fess 5 
525 1.1834 1.3764 

3 235 f.O778 0.7605 FX = -3657 
~ 45 1.0569 0.9966 Po —= 11069 
795 1.0090 1225972 

4 «35 1.3465 0.7216 FX = -4454 
o 45 Ieee ales 0.9272 FZ = 10666 
2.5 2207 SD i375 

5 oS 0.7109 Oi5 6 73 FX = -2103 
~ 45 Ombi/Z 3 Oy) 335 FZ = 10540 
255 0.6100 0.9169 

6 =25 1.3789 1.2997 FX = —-5982 
245 1.3247 16083 FZ = 16742 
5 1.2453 2.0279 

2. Finite Element Analysis Results 


a. Load Case 1 

In looking at the 275 knot load without static 
aeroelastic twist, it was noted that the largest stress 
concentrations in the structure were located in the rib legs, 
with the lower legs experiencing approximately 15 ksi in 
tension along their upper flanges and the upper legs seeing 
about -14 ksi (compression) along the lower flanges. These 
stresses were evenly distributed, as may be seen by the values 
in Figures 22, 23 and 24 where the major principal axis stress 
(o;) contours are shown. 


The minor principal (o,,) stresses, 


Von Mises and shear stresses show a similarly even 


distribution, with all stress levels well below the yield 
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stress value (Clemn 40 ksi) for the material. It should be 
noted that the apparent deformations in the plot are 
exaggerated and not to scale. This scaling provides the 
viewer with a better perception of the direction of 
displacement occurring, although in reality, the displacements 
are only on the order of 0.06 inches for this load case as 


determined by MSC/pal 2” 
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Figure 22. Major principal axis stress contours on 
inboard (WS256) end rib in untwisted 
condition (Case 1). 
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Figure 23. Major principal axis stress contours on 
middle (WS282) rib in untwisted condition 


(Case 1). 
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Figure 24. Major principal axis stress contours on 
outboard (WS320) rib in untwisted condition 
(Case 1). 
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The untwisted load case (case 1) does not 
accurately reflect the behavior of the leading edge since it 
does not include the effect of aeroelastic-induced spanwise 
twisting (torsion) of the wing. The stress contours in 
Figures 22, 23 and 24 were presented for comparison to load 
case 2d. 

b. Load Case 2 

A dramatic difference in the observed stress 
contours occurred when the finite element model was subjected 
to the net twist occurring in the wing box as determined by 
the static aeroelastic span load analysis. It should be noted 
that the displacement input to the model equated to only about 
0.3 degrees of front spar rotation from WS 256 to WS 320. The 
same set of stress contour plots as above are given in Figures 
25, 26 and 27. A striking contrast existed between the first 
and second load cases, with the case 2 major principal axis 
stress distribution very unevenly spread between the inboard 
and outboard ends. As may be seen in Figure 25, the inboard 
end rib (WS 256) experienced more than double the stress 
concentration in its lower leg with 34.8 ksi being the highest 
level contour shown. Moving outboard, WS 282 (Figure 26) 
showed a maximum of about 20 ksi in its lower leg while the 
stress in the same leg on the outboard (WS 320) rib went to 
zero. (A stress level of 2.9 ksi in the saddle and upper leg 


of the rib may still be seen.) 
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Figure 25. Major principal axis stress contours on 
inboard (WS256) end rib with twist applied 
(Case 2). 
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Figure 26. Major principal *axis stress contours on 
middle (WS282) rib with twist applied 
(Case 2). 
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Figure 27. Major principal axis stress contours on 


outboard (WS320) rib with twist applied 
(Case 2). 


Since it was observed that this same pattern of 
stress distribution occurred for the minor principal axis, Von 
Mises and maximum shear stresses, only the inboard end (WS 
256) contours are shown in Figures 28, 29 and 30. This 
pattern is to be expected since they are all geometrically 
related. The major and minor principal stresses occur on 
planes on which there is no shear stress and are oriented 
perpendicularly to each other, while the maximum shear stress 
occurs on planes which are at angles of 45 degrees to the 
principal planes. Von Mises stresses (o,) are derived from a 
criterion known as the Maximum Distortion Energy Criterion and 


may be found from the equation, 


magn uses 2 
Oy = OF O7r9r, + Oy 
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which is based on the determination of the energy associated 
with changes in shape of a given material. These relation- 
ships are valid under the assumption of a plane stress 
condition in the material. This means that the metal is thin 
in comparison to its other dimensions so that the stresses 
across the thickness of the metal may be considered as 
negligible. This plane stress condition is the case for most 
aircraft structural components since they are made of sheet 


material, and is valid here. [Ref. 19, 20] 
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Figure 28. Minor principal axis stress contours on 
inboard (WS256) end rib with twist applied 
(Case 2). 
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Figure 29. Von Mises Criterion stress contours on 


inboard (WS256) end rib with twist applied 


(Case 2). 
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Figure 30. Maximum shear stress contours on inboard 
(WS256) end rib with twist applied (Case 2). 





In Figures 31 through 36, the x and z displacements of 


the ribs are shown in pairs from inboard to outboard. 
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Note 


the small displacements occurring in the upper flange and 


lower hinge due to twist of the front spar. 
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Figure 32. Displacement (z) at inboard rib (WS256). 


67 


TRANS. DEFL. & 


A-3. 68B9E-82 
B-2. 7667E-82 
- 8444E-82 
. 2222E-83 
BBO0GE +48 
2222E-83 
8444E-82 
7667E-82 
68B9E-82 
6111 E-82 


bAUNrUawr 


C 
D 
E 
F 
G 
H 
I 
J 


TRANS. DEFL. 2 


. B66GE +48 
- @CLIE-82 
. 943B8E-82 
- JLSVE-82 
. 8B876E-82 
. B595E-82 
- 8631 E-@1 
- 2ABJE-B1 
- 41 75E-@1 
.9947E-B1 


ana = SECA we D 
pm pe ee ee 8 OT CD  @ 





Figure 34. Displacement (z) at middle rib (WS282). 
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Figure 35. Displacement (x) at outboard rib (WS320). 
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Figure 36. Displacement (z) at outboard rib (WS320). 


Another available representation of the stress 
concentrations present in the structure is an X-Y plot of all 
four stresses together in a form resembling a frequency 


scatter on an oscilloscope. This plot for the case 2 load 
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condition is presented in Figure 37. While this is a rather 
cluttered plot it does serve to provide, at a glance, the 
maximum and minimum stress concentrations in the structure. 
Since the stress contours of individual members in the 
structure have been examined and it has been found that the 
highest stresses are in the rib legs, the peak stresses 
depicted may be attributed to these locations. The positive 
(upper) portion of the plot is a combination of the major 
principal, Von Mises and shear stresses while the lower half 


depicts the minor principal axis stresses. For the sake of 
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Figure 37. X-Y plot of combined stresses (Case 2). 


brevity, and since the contour plots of the individual ribs 


retain the same relative form as those presented for cases 1 
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and 2, the stress levels for the remaining load cases will be 
presented in this form. 
c. Load Case 3 

For the 275-knot, 2.4-G rolling pullup, the 

maximum stress levels observed are approximately 33 ksi on the 

descending wing as shown in Figure 38. The loads on the 

leading edge of the ascending wing are lower than those for 

the descending wing. This speed was chosen because of the 

relationship which was shown to exist between leading edge 


loading, angle of attack and dynamic pressure. The result of 
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Figure 38. X-Y plot of combined stresses (Case 3). 
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higher airspeed at the same aircraft load factor may be see in 


part e of this Chapter. 
d. Load Case 4 
The 110,000 pound weight condition was chosen for 
this load case in order to examine the reduction in leading 
edge stress when operating within the normal flight envelope 
at less than maximum gross weight. The 240-knot airspeed is 
the approximate "corner" speed at a load factor of 3G for this 
weight. The maximum stress values seen under this condition 


are 28.8 ksi in tension and 28.3 ksi in compression as seen in 


Figure 39. 
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Figure 39. X-Y plot of combined stresses (Case 4). 
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e. Load Case 5 
In load case 5, the effect of a 2.4-G rolling 
pullup at 135,000 pounds gross weight at a true airspeed of 
350 knots was examined. In Figure 40, it may be seen that the 
maximum stress values are approximately 26 ksi, which is 
approximately 7 ksi less than the same maneuver at 275 knots 


(load case 3). 
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Figure 40. X-Y plot of combined stresses (Case 5). 


f. Load Case 6 
In order to examine the design ultimate flight 
load condition, load case 6 was run at 325 knots and 4.5G. 
The airspeed corresponds to the approximate stall speed 


(corner speed) for the 135,000 pound airplane gross weight. 


73 


Stress values of nearly 52 ksi may be seen in Figure 41. Note 
that this result exceeds the yield stress for the material (47 


ksi1). 
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Figure 41. X-Y plot of combined stresses (Case 6). 


3. Discussion 

The primary factor of concern in the observed loads . 
presented above is the apparently high level of stress in the 
rib legs as determined by the methods stated in the course of 
this writing. According to the Engineering Investigation (ETI) 
report for the USN P-3 mishap [Ref. 21] and a preliminary 
report issued by Aircraft Research Laboratories of Melbourne, 
Australia [Ref. 22], the initial structural failure in both 


mishaps was believed to have occurred in the lower rib legs, 
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with initial failure occurring in the outboard leg and 
proceeding inboard in a series of sequential failures. It 
would appear from the results of this study that this failure 
sequence may have been reversed, since the highest observed 
stresses were in the inboard rib. Additionally, the fact that 
there is an intermediate rib at WS 317, only three inches from 
the outboard end rib (rather than the eight- to nine-inch 
Spacing for all other ribs in the leading edge) supports the 
idea that the failure would likely have initiated at some 
other location. This close proximity of two ribs means that 
the load in that region would be shared, thereby reducing the 
stress in the outboard rib. In any case, the primary question 
remains as to whether there were sufficiently high stress 
levels in the rib legs to cause structural failure. 

While the maximum observed stress levels (found here 
for operation within the flight envelope, occurring under load 
case 2 in a wings level, 3-G pullup at 135,000 pounds gross 
weight) were some 12 ksi below the yield stress for the 
material, two areas of concern must be addressed: 1) the 
additional effect of extending this data from the assumed 64- 
inch model size to the true 92-inch leading edge segment which 
is installed on the aircraft; and 2) the effect of stress 


concentrations around holes in the rib legs. 


(Le 


a. Extending the Structure 

Proceeding under the assumption that the stress 
concentrations vary linearly with the length of the structure, 
extending the 64-inch model to the full (92-inch) length of 
the true structure would result in an increase in the observed 
stress in the lower rib leg to 51.9 ksi, which is in excess of 
the 47 ksi yield stress. Examining the major principal axis 
stresses in the lower legs of WS 256 and WS 282, this 
linearity would appear to hold true. Since aluminum is a 
ductile material, the stresses may subside as the material 
begins to yield, therefore, the structure may still carry the 
load without fracture. However, changes in the pressure 
distribution around the leading edge must also be considered 
as the structure begins to deform. Given the observed 
directions of the deformations as seen in the contour plots, 
it seems likely that the aerodynamic load would continue to 
increase due to the increased camber in the leading edge as it 
lifts perpendicular to the chord and develops a more circular 
form. This deformation, although initially small could add to 
the load, initiating a divergent stress scenario leading to 
failure. 

b. Stress Concentrations 

It is an established fact that stress 

concentrations around circular holes can multiply the 


localized stress by factors of between two and four, depending 
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on the directional orientation of the loading seen by the 
specimen. Along the inner surface of the end ribs is a 
continuous flange which is attached by a rivet spacing of 
approximately one-half inch. Additionally, there are four 
tapered angle stiffeners riveted to the web section of the 
rib. Two of these are symmetrically located on the upper and 
lower legs at a position approximately 52 percent of the 
leading edge chordwise dimension forward of the hinge. (On 
the outboard end rib (WS 320) this is approximately 12 inches 
forward of the hinge.) Around these rivet holes, stress 
concentrations may be assumed to occur. The occurrence of 
these concentrations at the onset of stress loading depends on 
the rivet installation procedure; §i.e., the stress 
concentration effect is delayed if the rivet hole is placed in 
compression (as in a dimpled or double dimpled installation 
feet 18)]). If the hole were initially in compression, a 
Margin or stress buffer would be provided, as the material is 
subjected to tension stress; i.e., the tensile stress applied 
must first exceed the compressive pre~stress of the rivet 
installation before the hole will begin to experience a stress 
concentration. The type of installation of the rivets under 
consideration is unknown by this author. Since the inside 
flange of the upper legs is placed under compressive loading 
nearly equivalent to that experienced by the lower leg in 
tension, the effect of compressive concentrations must also be 


considered. If the rivet holes were pre-stressed in 


fag) 


compression, this would add to the compressive stress level at 
this location, potentially causing the upper leg to fail first 
in compression. As before, this all depends on the stress 
reduction due to yielding which takes place in the ductile 
material. 

In the Australian incident report, the failure 
sequence was described as having initiated along the line of 
the riveted stiffener on the lower leg of the outboard end rib 
as mentioned above. In addition, it goes on to say that the 
failure then proceeded along a line of dimples in the lower 
legs of the intermediate ribs, from outboard to inboard. 
These dimples act as stiffeners for the intermediate ribs and 
would also act as stress concentrations under load. None of 
these stiffeners, on either outboard or intermediate ribs, nor 
their rivet holes, are included in the finite element model 
employed in this study. Again, the failure sequence described 
is opposite of that expected from the results of this study. 
The RAAF determination of the failure sequence was based on 
the observed outboard-to-inboard bending of the rib fragments 


which remained attached to the structure. [Ref. 22] 
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VI. CONCLUSIONS AND RECOMMENDATIONS 


A. CONCLUSIONS 

The results of this study indicate that failure of a P-3 
wing leading edge segment could be predicted to occur within 
the normal aircraft operating envelope providing proper 
account is given to the influences of wing static aeroelastic 
effects both upon wing span loads and torsional twists induced 
in the wing spar box. The location of the highest stresses 
are in the area of the observed failures as reported in the EI 
and Preliminary reports [Ref. 21 and 22]. It is reasonable to 
assume that stress concentrations around the rivet holes could 
be on the order of 1.5 or 2.0 times the observed stresses. 
This, combined with the effect of extending the model to the 
full 92-inch length, could result in stress levels in excess 
of the ultimate strength of the material, 60 ksi [Ref. 18]. 

The primary evidentiary conflict with this study is the 
reported location of apparent failure initiation at the 
outboard vice inboard end. The observed inboard bending of 
the rib fragments may have been induced by the upper portion 
departing the wing in some fashion other than that assumed. 
It is difficult, at best, to assess the direction in which the 


ribs would bend as the failure progressed. 
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While fatigue is potentially a contributor to this mode of 
failure, given the possibility of plastic strain having 
occurred earlier under a reduced load condition, the EI [Ref. 
21] stated that no undue hardness was detected in the failed 
USN structure. This would indicate that no material strain 
hardening had taken place up to the time of failure. 

Corrosion is another potential problem, particularly when 
dealing with aircraft operating routinely in low level 
overwater environments; however, no levels of corrosion which 
would effect the structural integrity of the leading edge were 
Clted in Ehe Er. 

Whether or not the aircrews overstressed the aircraft by 
exceeding the limitations of the flight envelope cannot be 
concluded here. Even they may not know for sure whether this 
was the case. Flight experiments conducted in the 2F87F P-3 
Simulator at NAS Moffett Field, CA, during the course of this 
analysis indicate that at high gross weights and aft center of 
gravity conditions, it is easy to exceed the 3G limit with 
just a firm pull on the control yoke. The rapid onset of G 
overload may be imperceptible due to the rate at which it was 
observed to occur. The cockpit G meter is not an instrument 


which is normally kept in the pilot’s scan. 
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B. RECOMMENDATIONS 

This analysis should be confirmed by independent 
validation. If these results are verified, re-enforcement of 
the end rib may be feasible; alternately, some limitations 
should be placed on the operational flight envelope of the 
P=3. Covered here are some recommended procedures for 
validation along with approximations of interim flight 
envelope limits, should validation and/or repair not be 
possible in the near-term. 

1. Validation 

Since the Naval Air Systems Command, through 

Aerostructures, Inc., has the base model from which this 
finite element model was derived, it is recommended that 
validation be conducted with the following modifications on 
the existing model: 

e Skin stiffness should be included in the model. The 
contribution of the skin to the load bearing properties of 
the structure is essential. Skin stiffness could be 
Simulated by artificially increasing the thickness to 
Simulate the double skin as was done here, or preferably, 
through addition of aw axis (quadrilateral element local 
coordinate system) stiffness factor which would emulate 
the combined stiffness of the inner and outer skins. 

e Static aeroelastic twist must be included in the analysis 
Since it was the factor contributing to the dramatic 
increase in inboard end rib stress. Since the NASTRAN” 
model is already of the appropriate length, this model 
would serve to verify or discount the theory that the 
stress concentrations would increase linearly with 


extended length. 


e A small finite element model of the area of immediate 
concern (the lower leg of an end rib) should be 


81 


constructed and analyzed. This model should be of 
sufficient detail to include the rivet holes and 
stiffeners in order to examine the amount of stress 
concentration occurring. Any pre-stress of the rivet 
holes should be taken into account. 

Laboratory tests of actual end ribs would prove 
helpful in determining the effect of the rivet holes on stress 
concentrations and provide data on the load levels required to 
cause fracture. 

2. Flight Envelope Modifications 

If these results are verified, the P-3 flight envelope 
should be restricted to prevent another mishap. If, in the 
estimation of NAVAIR, verification cannot be accomplished 
within a reasonable amount of time (six months), then it is 
recommended that an interim measure be taken to restrict the 
flight envelope until such verification or disproval can be 
accomplished. While the precise envelope restrictions would 
be developed by NAVAIR, it is recommended on the basis of the 
results observed here that the following limitations be 
applied if interim limitations are deemed appropriate: 

e Reduce the maximum sustained load factor for the aircraft 
from 3G to 2G at all operating weights above 100,000 
pounds. 

e Reduce the maximum sustained load factor for the aircraft 


from 3G to 2.5G at all operating weights up to and 
including 100,000 pounds. 
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While these are rough estimates, they are believed to provide 
a sufficient margin of safety to allow safe flight without 


severely impacting daily P-3 operations. 
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APPENDIX A 


ae Program "P3SPNLD3.BAS" 
Solve Static Aeroelastic Spanload Problem 
for a wing with straight taper 


Vortex Lattice capability uses swept bound vortex 


elements 
Inputseconsisteet: 
AR = Wing Aspect Ratio = B*2/S 
TR = Wing Taper Ratio = Ct/Cr 
SWP25 = 0.25 ch’d. sweep angle, +’ve is sweepback 
CEA = Elastic axis location on const. fraction 


chord line 
Subsonic Mach no. for aerodynamic compress. 
correctn. 
No. of equal length spanwise stations, RH wing 
No. of equal chordwise stations 
x* Comments ** 
Allows Vortex Lattice solution with MxN boxes on 
RH wing 
Max. of MxN = 50:’ Consistent with dimension 
statements 
N = 1 for elementary (Modified Weissinger) lifting 
line theory 
Wingspan, B, set to 1188 inches (for P3 applicatn. ) 
100 ’ SDYNAMIC 
DIM X1(50), Y¥1(50), X2(50), Y¥2(50), X3(50), Y3(50), 
SWP(50),DIM A1(50, 50), A2(50, 50), ASYM(50, 50), 
AANT(50, 50),DIM S(50, 50), F(20, 50), G(20, 50), 
XEA(20), YEA(20), EI(20), GJ(20),DIM SWPM(50), 
DELA(50, 20), MX(20, 50), MY(20, 50),DIM ALPHA(S0), 
U(20, 20), UGJ(20, 20), UEI(20, 20),DIM A(50, 50), 
XIN(50): ’ for Spanload solutns. using ELU/SLVB S/R’s 
DIM ALPHA1(50), ALPHA2(50), ALPHA3(50), ALPHA4(50), 
Q(50), DELCAM(S) 
OPEN "C:QBFILES\EXTRA.DAT" FOR OUTPUT AS #1 
150 ’ Input wing geometric information .. P3 usage 
AR = 7.539: TR = .40088: M = 10 
‘PRINT : PRINT "Aspect Ratio, AR ="; AR 
‘PRINT "Taper Ratio, TR ="; TR 
‘PRINT "No. R.H. Wing Spanwise Stas. ="; M 
"MACH = .6: PRINT TAB(5); "Mach No. ="; MACH 
INPUT "Airspeed (knots) ="; VEL 
MACH = (VEL * 1.68781) / 1116.3 


= 
lt at 


= = = = = = ™~ = = = ™= = = ~= = ™= ~ ™= ™= ™= ™= = ™ 
. 
Il 


PRINT "Mach No. =", MACH 
IF MACH > .95 THEN GOTO 151 
FC = 1! / SQR(1! ~- MACH * 2): GOTO 152 
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Ea loec = 1! 
152 ’ Continue Dummy statement 
N= 5 
‘PRINT "No. Chordwise Stas. ="; N 
‘INPUT "No. Chordwise Stas. ="; N 
SWP25 = -1.312: ‘’PRINT "C/4 Sweep (deg.) ="; SWP25 
‘INPUT "0.25 chd. sweep (deg.) ="; SWP25 
CEA = .4: ’PRINT "Elast. Axis Locatn., X/C ="; CEA 
‘INPUT "Elastic Axis Locatn., X/C ="; CEA 
‘INPUT "Dynamic Pressure, Q (psi) ="; Q 
INPUT "Is this a RIGID run? 1=Yes, 2=No, ANS=", ANS 
IF (ANS = 1) THEN 
® - 0 
ELSE 
Q = (.7 * 2116.2 * MACH * 2) / 144 
END IF 
PRINT "Q =", Q 
KMAX = M * N 
f 
‘ Calculate propeller thrust per engine (using curve fit 
‘ from 2F87F P-3 simulator data) 
THRUST = (5849 - 5.069 * VEL) 
‘ Calculate velocity boost due to prop slipstream 
v = VEL * 1.68781 
V1 = (-v + SOR((v * 2 + (4 * .05 * THRUST / (143.14 * 
.0023769))))) * 10 
‘ Use velocity boost to find boosted Q (Q1) 
Gu (85 * .0023769 * (v V1) * 2) / 144 
IF (ANS = 1) THEN 
Ol = 0 
END IF 
f 


‘ Set up Q as a vector with boosted pressure (Q1) in 
‘ region behind propellers 

FOR I = 1 TO KMAX 

Q(I) = Q: NEXT I 

“FOR I = 11 TO 35 

Q(I) = Ql: NEXT I 


PI = 4! * ATN(1!): ’ Establish constant 
SWP25 = SWP25 * PI / 180!: ’ Convert sweep to radians 


180 ’ Print header 


“GOSUB 1100 
"FOR MM = 1 TO 4 


200 ’ Determine Initial Wing Geometry 
| B= 1188!: ’ DEFAULT Value of P3 Wing Span, inches 
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DELB = .5 * B / M: ‘Spanwise spacing increments of 

: vortex elements 

CR = 2! * B / (AR * (1! + TR)): * Root Chord, anes 

S =B%* 2 / AR: ’ Wing Area, sq. in. | 

TANLE = TAN(SWP25) + (.5 * CR * (1! - TR)) / B: 

’ Tangent L.E. sweep 

'‘ Develop Mean Aero. Chord information 

MAC = 2! * CR * ((1! + TR) — (TR / (1! + TR))) / 32: 

‘ Cmac, inches 

IF TR = 1! THEN GOTO 210 

YMAC = .5 * B * (1! -— (MAC / CR)) / (1! - TR): GOTO 211 
210 YMAC = .5 * B 
211 XMAC25 = (YMAC * TANLE) + .25 * MAC 


250 ’ Determine coords. for wing vortex lattice corners & 
‘ €ontrol prc. 


CONST. = 1! = 3TR 

FOR I = 1TOM 

Cl = CR * (1! - (CONST1 * (I - 1!) / M)): ‘Total wing 

‘ ehord, imbd~svomcex 

C2 = CR * (1! - (CONSTI * IT /.M)):. % Wing chord, gene 0-8 
' vortex 

C3 =.CR.* (1! - (CONST1 * (L.-..5) /£.M)): * Choma 


‘ control point sta. 

DELC1 = C1 / NeeDELG2Q0= C27 Ne DELC 6 — con a 

YEA(I) = (I - .5) * DELB: ‘ Create "M" Elastic axis 
‘ coordinates 

XEA(I) = YEA(I) * TANLE + CEA * C3 

FOR J = 1TON 

K = (I~ 1) * N+ J: 'Create Vortex Lattice numbering 


‘ scheme 

Y1(K) = (I = 1) * DELB: Y2(h) = Yi@K) +SDEEE 
Y3(K) = Yi(K) + (.5 * DELB) 

X1(K) = DELC1 * (J - .75) + Y1(K) * TANLE 
X2(K) = DELC2 * (J - .75) + Y2(K) * TANLE 
X3(K) = DELC3 * (J - .25) + Y3(K) * TANLE 
TANSWP (X2(K) -~ XITCK)) / (Y2CK) =S¥ eee 


SWP(K) = ATN(TANSWP) 


TANSWPM = FC * TANSWP 
SWPM(K) = ATN(TANSWPM) 
NEXT J 
NEXT I 


‘ Tangent of Elastic Axis Sweep, etc. 

DELEA = (XEA(2) ~ XEA(1)): REA = SQR(DELEA ~ 2 + DELB“2) 
TANEA = DELEA / DELB: SINEA = DELEA / REA 
COSEA=DELB/REA 

WRITE #1, Q 

INPUT "ENTER G LOAD, n =", NZ 

‘Calculate required CL 

CLREQD = (NZ * 135000) / (.7 * 2116.2 * MACH * 2 * 1300) 
PRINT "CLREQD =", CLREQD 
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270 


300 


310 


Sal 


315 


WRITE #1, NZ 

‘ Bring in EI(I) and GJ(I) values for Structural model 
‘ I = 1 (Root Sta.) & I =M (Tip Sta.), M = 10 by 

‘ default 

‘ EI Data, Est. for P3 wing 

DATA 8.30E+10, 6.80E+10, 5.27E+10, 4.10E+10, 2.80E+10 
DATA 2.00E+10, 1.33E+10, 0.85E+10, 0.55E+10, 0.35E+10 
FOR I = 1 TO M: READ EI(I): NEXT I 

‘ GJ Data, Est. for P3 wing 

fee 7 SORTIO, 5.25810, 3.90E+10, 2.90E+10, 2.10E+10 
beeen 1,.30E+10, 0.85E+10, 0.65E1+10, 0.35Et10, 0.25E+10 
FOR I = 1 TO M: READ GJ(I): NEXT I 


‘Develop Aerodynamic Influence Coefficients 

‘ A1l(K1,K2) = RH Wing 

'’ A2(K1,K2) = LH Wing 

FOR K1 = 1 TO KMAX: ’ Kl Sta. is at the Control Point 

FOR K2 = 1 TO KMAX: ’ K2 Sta. is at the Vortex Station 
NUM1 = X3(K1) - X1(K2): NUM2 = X3(K1) - X2(K2) 

NUM1 = FC * NUM1: NUM2 = FC * NUM2 

‘ Bring in Prandtl-Glauert factor 

DEN1 = Y3(K1) - Y1(K2): DEN2 = Y3(K1) - Y2(K2) 

Rl = SOR(NUM1 ~*~ 2 + DEN1 * 2): R2 = SQR(NUM2*2+DEN2%2) 
‘ Find trig. functns for orthogonal transformation on 

‘ swept bound vortex 

SINSWP = SIN(SWPM(K2)): COSSWP = COS(SWPM(K2) ) 

H = NUM1 * COSSWP - DENI * SINSWP 

Y1ROT = NUM1 * SINSWP + DENI * COSSWP 

Y2Z2ROT = NUM2 * SINSWP + DEN2 * COSSWP 

COSTHETI1 = Y1ROT / R1: COSTHET2 = -Y2ROT / R2 

‘ Logic Check to avoid division by zero 

IF (ABS(H)) <= .001 THEN GOTO 310 

DELWBD = (COSTHET1 + COSTHET2) / H: GOTO 311 


DELWBD = 0O!: ‘ No downwash contributn. from bound 
i vortex 
‘’ Dummy statement space 


COS1 = NUM1 / R1: COS2 = NUM2 / R2 

DELWLH = (1! + COS1) / DEN1 

DELWRH = (1! + COS2) / DEN2 

Al(K1, K2) = (DELWLH + DELWBD - DELWRH) / (8! * PI) 
NEXT K2 

NEXT Kl 


‘ ** Similar logic for LH Wing Panel Aero. Influence 
: coeffs. 


FOR K1 = 1 TO KMAX: ’ Kil Sta. is at the Control Point 
FOR K2 = 1 TO KMAX: ’ K2 Sta. is at the Vortex Station 
NUM2 = X3(K1) - X1(K2): NUM1 = X3(K1) - X2(K2) 

NUM2 = FC * NUM2: NUM1 = FC * NUM1 

’ Bring in Prandtl-Glauert factor 
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DEN2 = Y3(K1) + Y1(K2): DEN1 = Y3(K1) + Y2(K2) 

R1 = SQR(NUM1 * 2 + DEN1 * 2): R2 = SQR(NUM2*2+DEN2“~2) 
‘ Find trig. functns for orthogonal transformation on 
‘ swept bound vortex 

SINSWP = -SIN(SWPM(K2)): COSSWP = COS(SWPM(K2) ) 

H = NUM1 * COSSWP - DENI * SINSWP 

Y1ROT = NUM1 * SINSWP + DENI * COSSWP 

YZROT = NUM2 * SINSWP + DEN2 * COSSWP 

COSTHET1 = YLROT / R13 9@0STHET2Z5—eeiZROe 7FRZ 

‘ Logic Check to avoid division by zero 

IF (ABS(H)) <= .001 THEN GOTO 320 


DELWBD = (COSTHET1 + COSTHET2) / H: GOTO 321 
320 DELWBD = 0O!: ‘ No downwash contributn. from bound 
. vortex 


321 ° Dummy statement space 
COS1 = NUM1 / R1: COS2 = NUM2 / R2 
DELWLH = (1! + COS1) / DENI 
DELWRH = (1! + COS2) / DEN2 
A2(K1, K2) = (DELWLH + DELWBD - DELWRH) / (8! * PI) 
ASYM(K1, K2) AIT(CK1, K2 wer eAZ@GKL KZ) 
AANT(K1, K2) A1L(KI, K2) @°A2Z( Ki, K2) 
NEXT K2 
NEAT Kl 


400 ’ Determine Struct. Infl. Matrix, based on test of Q > 0 
‘ Note: Q = 0 psi case is rigid wing situation 
IF Q > O! THEN GOSUB 3000 


‘ Skip the option for Sym. or Anti-Symm. Soln. 
‘GOTO 800: ‘’Symm. solution branch 


700 ’ Select Symmetric or Antisymmetric Solution 
INPUT "Symm. or AntiSymm. Prob., S/A"; P§$ 


710 IF P$ = "S" THEN GOTO 800 
IF PS = "A" THEN GOTO- 900 
GOTO 1000 


800 ’ Find Additional Loading 
‘ GOTO 805: ’ Branch to find effect of wing washout 
‘ Introduce Additional type of alpha 
INPUT "ENTER AOA IN RADIANS, AOA =", AOA 
FOR K = 1 TO KMAX: ALPHA(K) = AOA: NEXT K 
‘GOTO 805 
‘PRINT "Elastic wing, Alpha = ", ALPHA(1) 
WRITE #1, ALPHA(1) 


‘ Introduce alpha due to built-in geometric twist 
‘ 0.0 deg. at root, linearly to -2.5 deg. at tip 


WASH = -2.5 * PI / 180! 
FOR I = 1TOM 
ETA = (I - .5) / M: DELWASH = ETA * WASH 


FOR J = 1 TO N- 
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K = (I - 1) * N + J: ‘’ Create panel numbering 
ALPHA1(K) = DELWASH: NEXT J 

NEXT I 

‘GOTO 805 

I 

‘ Introduce alpha (rad.) due to eng./prop assbly. 
’ dead-wgt. 

DATA —0.0012, -0.0042, -0.0084, -0.0138, —-0.0189 
DATA -0.0240, -0.0318, -0.0366, -0.0366, -0.0366 
' Read dead-wgt. alphas one at a time.. 

FOR I = 1TOM 

READ DELALPH 

FOR J = 1TON 

K = (I -1)* N+ J: ’ Create panel numbering 
ALPHA2(K) = NZ * DELALPH: NEXT J 

NEXT I 

a 

’ Introduce alpha (rad.) due to camber 

DATA -.0525, -.0240, 0.0, .0420, .0730 

‘ Read in alphas due to camber at the root (WS = 0) 
FOR I = 1TON 

READ DELCAM(I) 

NEXT I 

‘ Calculate the fifty slopes (one per control point) as 
‘ airfoil varies from 14% to 12% thickness from root to 
Rep 

FOR I 1 TOM 

FOR J 1 TO N 

ro—o(r — 1je *eN + J 

ALPHA3(K) = DELCAM(J) * (1! - (I - .5) / (7 * M)) 
NEAT J 

NEXT I 

I 


‘ Introduce alpha (rad.) due to aileron float angle 
‘’ Calculate aileron float angle (rad) based on velocity 


' (kts) 

FLAIL = -(1.4548 - .0090025 * VEL + .0000444 * VEL * 2) 
* .01745 

FOR I = 1 TO KMAX 


ALPHA4(K) = O!: NEXT I 


ALPHA4(39) = .2 * FLAIL 

ALPHA4(44) = .2 * FLAIL: ALPHA4(49) = .2 * FLAIL 
ALPHA4(35) = .35 * FLAIL 

ALPHA4(40) = FLAIL: ALPHA4(45) = FLAIL: ALPHA4(50)=FLAIL 


‘ Combine angles by choice for net alpha distrib. 

INPUT "Enter type of run desired: O = Additional only, 
1= Washout only, 2 = Dead wt only, 3 = Camber only, 

4= Aileron float only, 5 = Complete solution, TYPE =", T 
WRITE #1, T 
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FOR I 1 TO KMAX 

IF (T = 0) THEN 

ALPHA(I) = ALPHA(I): ‘Additional loading selection 
ELSEIF (T = 1) THEN 

ALPHA(I) = ALPHA1(I): ’ Washout alpha selection 
ELSEIF (T = 2) THEN 

ALPHA(I) = ALPHA2(I): ’ Dead-wgt. induced alpha 

’ selection 


Il 


ELSEIF (T = 3) THEN 

ALPHA(I) = ALPHA3(I): ’ Camber alpha section 

ELSEIF (T = 4) THEN 

ALPHA(I) = ALPHA4(I): ’ Aileron float selection 

ELSEIF (T = 5) THEN 

ALPHA(I) = ALPHA(I) + ALPHA1(I) + ALPHA2(I) + ALPHA3(I) 
+ ALPHA4(I): ’ Complete wing selection 

END IF 

XIN(I) = ALPHA(I) 

FOR J = 1 TO KMAX 

A(I, J) = ASYM(I, J) - S(7,°9).* OGlpeenNE ed 

NEXT I 

f 


‘ Find Spanload Soln. L/Q from XIN(KMAX) 

GOSUB 2000 

’ S/R returns XIN(I) as L/Q vector of length KMAX 

LIFT = O!: MOMENT = O!: FOR I = 1 TO KMAX 

LIFT = LIFT + XIN(I) 

MOMENT = MOMENT + XIN(I) * (XMAC25 ~— .5 * (X1(I) + 
X2(I))): NEXT I 

‘ Find Lift & Moment Coefficient 

Cl = Lipt © 2 eee 


Q 
= 
Il 


MOMENT * 2! * DELB / (S * MAC) 

NP = .25 (GM © sen) 

CLTOT = CLREQD - CL 

PRINT "CL =": CL; ", CM ="* CM: ",CLTOT =": CLTOT. 
“PRINT “Neut. Pt. at ( tecmac )'") Nr tee hae 

GOSUB 1100: ’ Print Header 

CAVE = .5 * CR * (1! + TR) 

FOR I = 1 TOM: ’ Sum up L/Q at eta value to get CLC 
ETA = (I - .5) / M: Ghe = 0! 

FOR J = 1TON: K= (I - 1) * NT J 

CCUG —IGLG. + XIN(K) 

NEAT J 

‘ Find struct. twist due to airloads at front panel for 
a sta. nyt 


TWIST = O!: J2 = (I-1)*Nt1 
FOR K = 1 TO KMAX 
TWIST = TWIST + S(J2, K) * Q(I) * XIN(K): ’ Struct twist 


‘ due to airload 
NEXT K 


a2 


€3 = CR * (1! — “CONSTI * (1 — .5) / M)): ‘’ Chord at 
eenerol point sta. 

SiseCr = CLC / C3 

erer="CECe/ (CL Se CAVE) 


GOSUB 1210: ’ Print eta, cLc/CLCave, cL(sect) & struct. 
twist 

NEXT I 

GOTO 990 


900 ’ Anti-Symmetric Solution Branch 
’ Set up alpha for roll damping due to pB/2V = HELIX 
INPUT "Enter run desired: O=Roll helix, 


1=Aileron deflection, RUN=", R 
’ Calculate available Pb/2V (rad) at specified velocity 
’ (kts) 
HELIX = (.21009 - .0008744 * VEL + .000001 * VEL ~* 2) 
IF (R 1) GOTO 905 


FOR I = 1TOM 

ALPDAMP = (I - .5) /M 

’ Multiply alpha for roll damping by available roll 
’ helix angle (rad) at given velocity (kts) 

ALPDAMP = ALPDAMP * HELIX 

FOR J = 1TON 

K = (I - 1) * N + J: XIN(K) = ALPDAMP: NEXT J 

NEXT I 

GOTO 908 


905 ’ Use following statement for aileron control 
’ effectiveness 
FOR I = 1 TO KMAX: XIN(I) = Of: NEXT I 
’ Calculate aileron deflection (rad) available at given 
’ velocity (kts) 
Deiat L = -(90.30699 — .65754 * VEL + .001771 * VEL *~ 2 - 
PeOWO0 16 *PVEL ~ 3) .01745 
’ Deduct float angle from available control wheel 
’ aileron deflection 


FLAIL = -(1.4548 —- .0090025 * VEL + .0000444 * VEL ~* 2) 
* .01745 

DELAIL = DELAIL - FLAIL 

MXIN(39) = .2 * DELAIL: XIN(44) = .2 * DELAIL 

XIN(49) = .2 * DELAIL 

XIN(35) = .35 * DELAIL 

XIN(40) = DELAIL: XIN(45) = DELAIL: XIN(50) = DELAIL 


908 ’ Solve Spanload Problem 
FOR I = 1 TO KMAX: FOR J = 1 TO KMAX 


A(I, J) = AANT(I, J) - S(I, J) * QO(1): NEXT J 
NEXT I 

’ Use S/R ELU/SLVB to solve anti-symm L/Q 
GOSUB 2000 


ROLL = 0!: FOR I = 1 TO KMAX 


a6 


ROLL = ROLL + XIN(I) * Y¥3(I): NEXT I 

CROLL = -ROLL * 2! * DELB / (S * B) 

IF (R = 0) THEN 

PRINT "Roll Damping Deriv. CLP ="; CROLL / HELIX 

ELSE 

PRINT "Aileron Control Deriv., CLa =", CROLL / DELAIL 
END IF 


910 ‘IF Q > 0 THEN GOTO 915 
‘CROLLQO = CROLL 


915 ’ Dummy Skip Statement 
‘EROLL = CROLL / CROLLQO 
‘GOSUB 1200: ’ Print results 
‘ Branch for cLc,CLsect,Twist 
GOSUB 1100: ’ Print Header 
CAVE = .5 * CR * (1! + TR) 
FOR I =1TOM ‘’Sum up L/Q at ETA to get CLC 
ETA = (I - .5) / M: CLC = 0! 
FOR J =1TO0N: K=(I-1) *N+JI 


CLC = CLC + XIN(K) 

NEAT J ' 

’ Find struct. twist due to airloads at front panel for 
a sta. my i 

TWIST = O!: J2 =(I-1)* NH+1 

FOR K = 1 TO KMAX 

TWIST = TWIST + S(J2, K) * Q(I) * XIN(K) 


‘ Struct twist due to airload 

NEAT K 

C3 = CR * (1! - (CONST1 * “I - 35) /SD ce Canora 

‘ control point sta. 

CUSECT ele. 7 63 

‘CLC = CLOWMCCL * (CAVE) 

GOSUB 1210: ’ Print eta, cLc, cL(sect) & struct. twist 


NEXT I 

‘GOSUB 1100 

"CAVE. = .5 4eCR * Gi! ear 

‘FOR I = 1TOM: ’ Defaulted for N = 1 (Chrdwse sta.) 
‘ETA = (I - .5) / M 

‘CLC = XIN(I) / CAVE: GOSUB 1200 


"NEXT I 


990 "ee —= Crt 1) 
"NEXT MM 
‘GOTO 700 


1000 END 
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1200 


1210 


2000 


2050 


’ S/R to Print header, whole selection shown 

‘PRINT TAB(3); "Mach No."; TAB(14); "CLalpha"; TAB(25); 
, "CmCl"; TAB(35); "N.P."; TAB(45); "Clp" 
Sere TAB( 5S); "O, psi"; TAB(14); "CLalpha"; TAB(25); 

t "CmCl";> TAB(35); " N.P."; TAB(45); 


r "Rlift": TAB(55); "Del-NP" 

‘PRINT TAB(10); "Roll Control Evaluation, Mid Ail. at 
: O05! 

‘PRINT TAB(5); "Q, psi"; TAB(15); "Croll"; TAB(25); 

: MErol* 

PRINT 


Poel TAB(6); “2Y/B"; TAB(13); "Clc/CLcave"; TAB(26); 

: "Cl"; TAB(35); "Twist"; TAB(45); "Degrees" 

RETURN 

‘ S/R to Print results, whole selection shown 

‘PRINT USING "#####.####"; Q; CL; CM; NP; ELIFT; DELNP 

‘PRINT USING "#####.####"3; Q; CROLL; EROLL 

PRINT USING "#####.####"; ETA; CLC; CLSECT; TWIST; 
TWIST / .01745 

WRITE #1, ETA, CLC, CLSECT, TWIST, TWIST / .01745 

RETURN 


‘ S/R ELU 

‘ Tri-Diagonalizes the input matrix A(KMAX,KMAX) 

‘ Input Column vector, XIN(KMAX), gets replaced by 
‘ the output which is returned to the main program 
NM1 = KMAX - 1 

FOR K = 1 TO NM1: KP1l = K+41 

FOR I = KP1l TO KMAX 

he= A(T, K) / A(K, K): A(I, K) =G 

FOR J = KP1l TO KMAX 

Pye, J) = A(l, J) +.G * A(K, J): NEXT J 

NEXT I 

NEXT K 

GOSUB 2050: ’ Use S/R SLVB for next step 

RETURN 


' S/R SLVB 
‘ Solves the Tri-Diagonalized matrix [A] obtained 
‘ from ELU by back substitution 


NM1 = KMAX - 1: NP1 = KMAX + 1 
FOR K = 1 TO NM1 

KPl =K+1 

FOR I = KP1 TO KMAX 

XIN(I) = XIN(I) + A(I, K) * XIN(K) 
NEXT I 

NEXT K 


XIN(KMAX) = XIN(KMAX) / A(KMAX, KMAX) 
FOR K = 2 TO KMAX 

I = NP1 - K 

Ji =I+1 


2: 


FOR J = J1 TO KMAX 

XIN(I) = XIN(I) - A(I, J) * XIN(J) 
NEXT J 

XIN(I) = XIN(I) / A(I, I) 

NEXT K 

RETURN 


3000 ‘S/R to Determine [S] Structural Infl. Coeff. matrix 
Ff 


3010 ’ Generate [My] and [Mx] matrices, order M x KMAX 
‘ First get diagonal elements 
K = 1: FOR I = 1TOM: FOR J =1TON 
MY(I, K) = .5 * (XEA(I) -— X1(K) - .75 * DELB * 
TAN (SWP(K) )) 
MX(I, K) = .25 
K = K + 1: NEXT J 
NEXT I 
‘ Now get Upper triangular elements 
VALUE = 1!: FOR L = 1tTOM-1 
K=L* N+ 1: K2 =M-L 
FOR I = 1 TO K2: FOR J = 1TON 


MY(I, K) = XEA(I) - X1(K) - .5 * DELB * TAN(SWP(R)) 
MX(I, K) = VALUE: K = K + 1: NEXT J 
NEXT I 


VALUE = VALUE + 1!: NEXT L 

‘ Test circuit 

PRINT "Print out of MX(I,J)" 

FOR I = 1 TO M: FOR J = 1 TO KMAX 
PRINT USING "##.##"; MX(I, J); 
NEXT J 

PRINT : NEXT I 


~ ~ ~ ™= = 


3020 ’ Generate [DELA(I,J)], order KMAX x M 
K = 1 
FOR I = 1TO-M: FOR J = 1TON 
DELA(K, I) = 1! 
Koon aes NEXT J 
NEAT I 
’ Testaciceuit 
PRINT "Print out of DELA(I,J)" 
FOR I = 1 TO KMAX: FOR J = 1TOM 
PRINT USING "##.##"; DELA(I, J); 
NEAT J 
PRINT : NEXT I 


™~ ™~ ~ ~ ™~ 


3030 ’ Generate [U(I,J)] Square matrix, order M x M 


FOR I = 1TO M: FOR J =1TOM 
IF I = J THEN U(I, J) = .5 
IF I < J THEN U(I, J) = 0! 


IF I > J THEN U(I, J) 1! 


916 


NEXT J 

NEAT I 

‘ Test circuit 

FOR I = 1TOM: FOR J = 1TOM 
PRINT USING "##.##"; U(I, J); 
NEXT J 

PRINT : NEXT I 


™~ ~™ ™ ™ 


3040 ’ Generate UGJ(I,J) & UEI(I,J) matrices, order MxM 
’ Embed DELB”2 type of constants 
DELB2 = DELB ~ 2: TANDELB2 = TANEA * DELB2 
FOR I = 1TO M: FOR J = 1TOM 
‘ DELB2 = 1!: GJ(J) = 6! - J: ’ for test circuit 


Meg(l, Jia DELB2 * U(I, J) / GJ(J) 
UEI(I, J) = TANDELB2 * U(I, J) / EI(J): NEXT J 
NEXT I = TANDELB2 * U(I, J) / EI(J): NEXT J 


‘Test Circuit 

’ PRINT "Sample of UGJ(I,J)" 

’ FOR I = 1 TOM: FOR J = 1 TOM 
‘ PRINT USING "##.##"; UGJ(I, J); 
’ NEXT J 

“PRINT : NEXT I 


3050 ’ Start gathering terms for [S] matrix assembly 
DELSIN = DELB * SINEA: DELCOS = DELB * COSEA 
FOR I = 1 TO M: FOR J = 1 TO KMAX 


Hated) = O08 Gil, J) = 0! 
FOR K = 1 TOM 
Bil, J) = F(L, J) + UGJ(I, K) * (COSEA * MY(K, J) + 


DELSIN * MX(K, J)) 

G(I, J) = G(I, J) + UEI(I, K) * (SINEA * MY(K, J) - 
DELCOS * MX(K, J)) 

NEXT K 

NEXT J 

NEXT I 

3060 ’ Form the [S] matrix, order KMAX x KMAX 
FOR I = 1 TO KMAX: FOR J = 1 TO KMAX 


cei J) = Ot 

FOR K = 1 TOM 

elie )je-ee(l, J) + DELA(I, K) * (F(K, J) + G(K, J)): 
NEAT K 

NEXT J 

NEXT I 

‘ {S] matrix is now available 

RETURN 
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APPENDIX B 


P-3 LEADING EDGE FAILURE ANALYSIS 
TWO DIMENSIONAL PANEL METHOD 


PROGRAM PANEL 


THIS FORTRAN PROGRAM ITERATIVELY COMPUTES THE PRESSURE 
AND FORCE DISTRIBUTION AROUND THE P-3 AIRFOIL AT 24 
SELECTED WING STATIONS AND WRITES A LOAD FILE FOR USE 
IN MSC PAL2 FINITE ELEMENT ANALYSIS SOFTWARE. INPUT 
CONSISTS OF THE AIRFOIL COORDINATES (X,Y) FROM A FILE 
CALLED P3.DAT, THE 24 WING STATION LOCATIONS IN 

INCHES AND THE NODE NUMBERS FROM THE FINITE ELEMENT 
MODEL AT WHICH FORCES AND CONSTRAINTS WILL BE APPLIED. 
OUTPUT CONSISTS OF THE PAL2 LOAD FILE (LE.LD), A 
TABULAR FILE IDENTIFYING THE COEFFICIENT OF PRESSURE, 
DYNAMIC PRESSURE AND FORCES AT A GIVEN PANEL MID POINT 
ON THE LEADING EDGE OF THE AIRFOIL. SECTION LIFT 
COEFFICIENT (Cl) IS MATCHED TO AN INPUT Cl REQUIREMENT 
(GENERATED BY AN AEROELASTIC SPAN LOAD ANALYSIS 
PROGRAM) THROUGH ITERATION OF ANGLE OF ATTACK. 


PARAMETER (N=53,M=N+1, PI=3.14159265, RHO=. 0023769, 
PAMB=14.7) 

REAL X(M),Y(M),THETA(M),BETA(M,M),R(M,M) ,XM(M),YM(M), 
A(100,101),AN(M,M),BN(M,M),AT(M,M),BT(M,M), 
SUMBN(M,M) ,SUMB1,SUMB2,SUMA(M),SUMB(M) ,VTAN(M), 
Q(M),CP(M),P(M),NUM(M) ,DEN(M) , ALPHA,V,VEL,AOA, 
CY(M),CX(M),CM(M),FY(M) ,FX(M) , VTOT(N) , CYT, CXT, CMT 
FYT, FXT,CD,CL,CLRQD(100),WS(100),XML(M),YML(M), 
LOADM(M) , INTWIST, OUTWIST, TWIST, XTOPDIS, XBOTDIS, 
XDISTOP(48),XDISHIN(24),ZDISTOP(48) , ZDISHIN(24) 

INTEGER K,1I,J,NODE(100,100),NTOP(100),NHIN(100) 


OPEN (UNIT=51, FILE=’'P3.DAT’ , STATUS=’ UNKNOWN“ ) 
OPEN (UNIT=57, FILE=’DUMP. DAT’ , STATUS=’ UNKNOWN “ ) 
OPEN (UNIT=58, FILE=‘’WS.DAT’ , STATUS= ’ UNKNOWN’ ) 
OPEN (UNIT=60, FILE=’LE.LD’ , STATUS=’ UNKNOWN“ ) 

OPEN (UNIT=61, FILE=’TOP.DAT’ , STATUS=’ UNKNOWN’ ) 
OPEN (UNIT=62, FILE=’HINGE.DAT’ , STATUS=’ UNKNOWN “ ) 
OPEN (UNIT=63, FILE=’NODESa. DAT’ , STATUS=’ UNKNOWN’ ) 
OPEN (UNIT=64, FILE=’NODESDb. DAT’ , STATUS=’UNKNOWN ’’ ) 
OPEN (UNIT=65, FILE=’NODESc.DAT’ , STATUS=’ UNKNOWN’ ) 
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C....-READ IN AIRFOIL DATA 


DO 10 K=1,N 
READ(51,500) X(K),Y(K) 
10 CONTINUE 
500 FORMAT(2F10.6) 


X(N+1)=X(1) 


¥(N+1)=Y¥(1) 
C.....READ IN FINITE ELEMENT NODE NUMBERS FOR LATER USE 
Cc (READ AS NODE(PANEL NO. ,RIB) ) 


DO 12 I=17,36 
READ(63,575) (NODE(I,J),J=1,8) 
READ(64,575) (NODE(I,J),J=9,16) 
READ(65,575) (NODE(I,J),J=17,24) 
12 CONTINUE 
575 FORMAT(9I5) 


ae WRITE HEADER FOR FINITE ELEMENT LOAD FILE 
WRITE(60,590) 

Steals « REQUEST AND RECEIVE AIRSPEED INPUT 
PRINT*, ENTER VELOCITY IN KNOTS (XXX.XX) ’ 


READ* , VEL 
V = VEL*1.6878 


GA Oe REQUEST AND RECEIVE FRONT SPAR TWIST ANGLES 
PRINT*, ’ENTER TWIST AT 2Y/B=.35 (rad) (.XXXX) ’ 
READ* , TWIST35 
PRINT*, ’ENTER TWIST AT 2Y/B=.45 (rad) (.XXXX) ’ 
READ* , TWIST45 
PRINT*, ‘ENTER TWIST AT 2Y/B=.55 (rad) (.XXXX) ’ 
READ* , TWISTSS 
Cae. - CALCULATE CHORD LENGTH OF P-3 WING AT SELECTED WING 
C STATIONS 
Go oe SET UP OUTER LOOP TO ITERATIVELY CALCULATE PRESSURES 
c FOR A SERIES OF WING STATIONS 
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Or INPUT WING STATIONS, Cl AND LOAD MULTIPLIER FROM FILE 
C WS.DAT AND CALCULATE CHORD LENGTH AT THE WING STATION 


502 FORMAT(F7.1,F8.4,F6.2) 
DO 15 L=1,24 
READ(58,502) WS(L),CLROQD(L),LOADM(L) 
PRINT*, /ITERATION=’,L 
PRINT*, ‘WS=’,WS(L) 
PRINT*, ’CLRQD=’ , CLROD(L) 
CHORD = 227.*(1-(0.6*(2.*WS(L)/1188.))) 
PRINT*, ‘CHORD=’ , CHORD 


Ca CALCULATE THE THICKNESS FRACTION AT A GIVEN WING 
6 STATION (USE TO LINEARLY VARY AIRFOIL THICKNESS WITH 
C WING STATION) 
DTHICK = 1-((1-(12./14.))*(2.*WS(L)/1188.)) 
PRINT*, ‘THICKNESS RATIO =’ ,DTHICK 
Cee. . SCALE THE THICKNES OF THE AIRFOIL AT EACH WING STATION 
E BY THE THICKNESS FRACTION 


DO 16 K=1,N 
Y(K)=Y(K) *DTHICK 
16 CONTINUE 


Cay «ae SET STARTING ALPHA 
ALPHA = (PI/180. ) 
OF SET COUNTER START FOR AOA ITERATION 
KOUNT = 1 
C....e-BEGIN ITERATION OF ALPHA TO MATCH CLREQD = CL 


1100 CONTINUE 
omen. FIND NORMALIZED MID-POINT OF PANEL 
DO 30 I=1,N-1 
XM(I)=0.5*(X(1I)+X(I+1)) 
YM(I)=0.5*(¥(1I)+Y(I+1)) 
CC. .e FIND MID-POINT OF PANEL IN INCHES 


XML (I)=XM(I)*CHORD 
YML(1I)=YM(I) *CHORD 
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C.....-CALCULATE THETA (ANGLE BETWEEN PANEL AND CHORD LINE) 


NUM(I)=Y¥(I+1)-Y(I) 
DEN(I)=X(I+1)-X(I) 
THETA (I)=ATAN2(NUM(I),DEN(I)) 


C....-CALCULATE R (DISTANCE BETWEEN PANEL MID-POINTS) 


DO 20 J=1,N 
Perr SORT CCRT eas) per ii LY (0) ez) 


C....e-CALCULATE BETA (ANGLE AT MID POINT OF PANEL I, 
Cc ENCLOSING END POINTS OF PANEL J) 


BETA(I,J)=ATAN2(((YM(I)-Y(J+1))*(XM(I)-X(J)) 
# 1 1) er aed COIL (I) 
# ((XM(I)-X(T+1) ) *(XM(1T)-X( J) ) 
# patil) Y (i ele) CM (1) SY (oa) 


20 CONTINUE 
30 CONTINUE 


Oe CALCULATE AUGMENTED "A" MATRIX 


SUMB1=0.0 
SUMB2=0.0 
DO 50 I=1,N-1 
SUMBN(I,N)=0.0 
DO 40 J=1,N-1 
IF (I.EQ.J) THEN 
AN(I,J)=0.5 
AT(I,J7)=0.0 
BN(I,J)=0.0 
BT(I,J)=0.5 
ELSE 


AN(I,J)=(SIN(THETA(I)-THETA(J))*LOG(R(I,J+1)/R(I,3J)) 
# +COS (THETA(I)-THETA(J))*BETA(I,J))/(2*PI) 


BN(1I,J)=(COS(THETA(I)-THETA(J) )*LOG(R(I,J+1)/R(I,3J) ) 
# ~SIN(THETA(I)-THETA(J))*BETA(I,J))/(2*PI) 

AT(I,J)=-BN(I,J) 
BT(I,J)=AN(I,J) 

END IF 

SUMBN(I,N)=SUMBN(I,N) + BN(I,J) 

IF (I.EQ.1) THEN 
SUMB1=SUMB1 + BT(1,J) 

ELSEIF (I.EQ.N-1) THEN 
SUMB2=SUMB2 + BT(N-1,J) 

END IF 
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A(I,d) = ANCL oD 
A(I,N) = SUMBN(I,N) 
A(N,J) = AT(1,J) + AT(N-1,J) 
A(N,N) = SUMB1 + SUMB2 

40 CONTINUE 


50 CONTINUE 


DO 60 I=1,N-1 
A(I,N+1) = -V*SIN(ALPHA-THETA (LI) ) 
60 CONTINUE 


A(N,N+1)=—-V* (COS (ALPHA-THETA(1) )+ 
+ COS (ALPHA-THETA(N-1) ) ) 


CALL GAUSS(N,A) 


GAMMA=0 . 0 
DO 80 I=1,N 
QO(I)=A(1I,N+1) 
GAMMA=A(N,N+1) 
80 CONTINUE 


DO 100 I=1,N-1 
SUMA(I)=0.0 
SUMB(I)=0.0 
DO 90 J=1,N-1 
SUMA(I)=SUMA(I)+Q(J)*AT(I,J) 
SUMB(I)=SUMB(1)+BT(I,J) 
VTOT (I )=SUMA(I)+GAMMA*SUMB(I)+V*COS (ALPHA-THETA (I) ) 
VTAN(I)=VTOT(I)/V 
CP(I)=1.0-(VTAN(I))**2.0 
P(I)=((CP(I)*0.5*RHO*(V**2.))/144.) 
90 CONTINUE 
100 CONTINUE 


C.....CALCULATE FORCE COEFFICIENTS (CX,CY) AND FORCES (FX= 
Ee CHORDWISE, FY=NORMAL TO CHORD) AT PANEL MID POINT 


CYT=0 

CXT=0 

CMT=0 

FXT=0 

FYT=0 

DO 110 I=1,N-1 
CY(I)=-1.0*CP(I)*(X(I+1)-X(I)) 
CX( Wie CP (1) * Grr yr) 


C sjesmene MULTIPLY BY CHORD FOR FULL SCALE SOLUTION 


FY (I)=CHORD*CY(I)*0.5*RHO*(V**2.)/144 
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FX (I) =CHORD*CX (I) *0.5*RHO*(V**2.)/144 
Cr) Sele) Genie ACI) ) SAM rl) + 
+ Cael iat Lp yeyM¢ 1 )) 

CYT=CYT+CY (TI) 

CXT=CXT+CX (I) 

FYT=FYT+FY (I) 

FXT=FXT+FX (I) 

110 CONTINUE 


C.....GET FORCE TOTALS FORWARD OF FRONT SPAR (.15C) 


FYNT=0.0 
FXNT=0.0 
DO 115 I=18,35 
FYNT=FYNT+FY (TI) 
FXNT=FXNT+FX (TI) 
115 CONTINUE 


C....-CALCULATE APPROXIMATE TOTAL FORCES ON COMPLETE LE 
8 SECTION (BETWEEN NACELLES) BY TAKING FORCES AT 
C ree POINT * 92 INCHES 


IF (L.EQ.7) THEN 
FYTOT=FYNT*92 
FXTOT=FXNT* 92 

ENDIF 


Care! PERFORM AXIS ROTATION TO FIND CL 
CL=CYT*COS (ALPHA) -CXT*SIN( ALPHA) 
Cae... PERFORM ITERATION OF ANGLE OF ATTACK 


DIFF=.0001 

CLDIFF=CLROD(L)-CL 

KILL=50 

DELTA=0.1745 

IF (KOUNT.GT.KILL) THEN 
PRINT*, ‘COUNTER EXCEEDED’ 
GOTO 111 

ELSEIF (CLDIFF.GT.DIFF) THEN 
PRINT’ (I3,2F10.6)’,KOUNT, ALPHA, CL 
KOUNT=KOUNT+1 
ALPHA=ALPHA+ (DELTA*ABS(CLDIFF) ) 
GOTO 1100 

ELSEIF (CLDIFF.LT.0.0) THEN 
PRINT’ (13,2F10.6)’,KOUNT, ALPHA, CL 
KOUNT=KOUNT+1 
ALPHA=ALPHA- (DELTA*ABS(CLDIFF) ) 
GOTO 1100 

ELSE 


nO Ss 


PRINT’ (4F10.5)’,ALPHA, CL, FXNT, FYNT 
WRITE(57,570) VEL,WS(L),ALPHA, CLRQD(L),CL,FXNT, FYNT 
WRITE(57,560) 

ENDIF 


Cw cece WRITE PRESSURES AND FORCES TO OUTPUT FILES 


DO 120 I=17,36 
WRITE(57,580) I,XM(I),XML(I),YML(1I),CP(1),P(I), 
A FX (I) *LOADM(L) ,FY(1I)*LOADM(L) 
120 CONTINUE 


Cw cc ee WRITE LOADING FILE (LE.LD) FOR FINITE ELEMENT 
APPLICATION (PAL2). 


FX(17)=FX(17)*0.5 
FY (17)=FY(17)*0.5 
FX (23)=FX(23)+0.5*FX(24) 
FY (23)=FY(23)+0.5*FY(24) 
FX (25 )=FX(25)+0.5*FX(24) 
FY(25)=FY(25)+0.5*FY(24) 
FX (26 ).=FX (26 )+FX(27) 
FY (26)=FY(26)+FY(27) 
FX (28)=FX(28)+0.5*FX(29) 
FY(28)=FY(28)+0.5*FY(29) 
FX (30) =FX(30)+0.5*FX(29) 
FY (30)=FY(30)+0.5*FY(29) 
FX (36)=FX(36)*0.5 
FY (36)=FY(36)*0.5 
DO 125 I=17,36 
J=I 
IF (I.LE.23) THEN 
WRITE(60,591) NODE(J,L),FX(I)*LOADM(L), 
an NODE(J,L),FY(1I)*LOADM(L) 
ELSEIF (I.EQ.25) THEN 
WRITE(60,591) NODE(J,L),FX(1I)*LOADM(L), 
a NODE(J,L),FY(1)*LOADM(L) 
ELSEIF (I.EQ.26) THEN 
WRITE(60,591) NODE(J,L),FX(1I)*LOADM(L), 
+ NODE(J,L),FY(1I)*LOADM(L) 
ELSEIF (1I.EQ.28) THEN | 
WRITE(60,591) NODE(J,L),FX(1)*LOADM(L), 
+ NODE(J,L),FY(1)*LOADM(L) 
ELSEIF (1.GE.30) THEN 
WRITE(60,591) NODE(J,L),FX(I)*LOADM(L), 
+ NODE(J,L),FY(1)*LOADM(L) 
ENDIF 
125 CONTINUE 
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Soe WOnMAi ZX, “PANEIG’, 3X, “AMCE)’,4X, ‘AML(1I)”,4X, ‘YML(1)’, 


+ 5X, 'CP(I)’,5X,’P(1)’,6X, 'FX(1I)’,5X, 'FY(I)") 
570 FORMAT(//,2X, 'VEL=’,F10.1,/,2X, 'WS=',F10.3,/,2X, 

+ ‘ALPHA=’,F10.5,/,2X, 'CLRQD=’,F10.5,/,2X, 'CL=’, 

+ F10.5,/,2X, 'FXNT=',F10.5,/,2X, 'FYNT=',F10.5,/) 


580 FORMAT(2X,13,7F10.5) 
590 FORMAT(1X, 'FORCES AND MOMENTS APPLIED 0’) 
591 FORMAT(1X, ’FX’,17,F10.5,/,1X, 'FZ’,17,F10.5) 
111 CONTINUE 

15 CONTINUE 


C....--e-CALCULATE TWIST DISPLACEMENTS OF FINITE ELEMENT NODES 


INTWIST=0.82155*(TWIST45-TWIST35 )+TWIST35 
OUTWIST=0 .89371* (TWIST55-TWIST45 )+TWIST45 
TWIST=OUTWIST-INTWIST 

XTOPDIS=TWIST*8.8 

XBOTDIS=(-1.0)*XTOPDIS 


DO 130 I=1,24 
J=I-1 
XDISTOP(2*I-1)=(XTOPDIS/23.)*J 
XDISTOP(2*1I)=(XTOPDIS/23.)*J 
XDISHIN(1I)=(XBOTDIS/23)*J 
TT=(TWIST/23.)*J 
XT=(XTOPDIS/23.)*J 

IF (TT.LE.0625) THEN 
ZDISTOP(2*I)=0.5*TT*XT 

ELSE 
ZDISTOP(2*1I)=(-0.5)*(TT-.125)*XT 

END IF 
ZDISTOP(2*I-1)=(-0.5) *TT* (XT+1.0) 
ZDISHIN(I)=0.5*TT*XT 

130 CONTINUE 


WRITE(57,555) FXTOT,FYTOT | 
555 FORMAT(/,2X,’APPROX TOTAL FX(92")=’,F10.1, 


a /,2X,'APPROX TOTAL FY(92")=’,F10.1) 
oS ore READ AND WRITE NODES AND CORRESPONDING BOUNDARY 


c CONDITIONS IN FINITE ELEMENT LOAD 


WRITE (60,597) 


READ(61,592) (NTOP(I),I=1, 48) 
READ(62,601) (NHIN(I),I=1,24) 
WRITE(60,593) (NTOP(I),I=1,48) 
WRITE(60,5931) (NTOP(I),XDISTOP(I),I=1,48) 
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Bee ENFORCE Y-AXIS DISPLACEMENT CONSTRAINT ALONG FRONT 
EDGE OF SPAR CAP FLANGE ONLY. (RELEASE BACK EDGE TO 
SIMULATE SINGLE LINE OF SCREW FASTENERS. ) 


WRITE(60,5932) (NTOP(2*I),I=1,24) 
WRITE (60,5933) (NTOP(I),ZDISTOP(I),I=1,48) 
WRITE(60,594) (NHIN(I),XDISHIN(I),I=1,24) 
WRITE(60,599) (NHIN(I),I=1,24) 
WRITE(60,600) (NHIN(I),ZDISHIN(I),I=1,24) 
WRITE(60,595) (NHIN(I),I=1,24) 
WRITE(60,596) (NHIN(I),I=1,24) 
WRITE (60,598) 
592 FORMAT(IS) 
601 FORMAT(IS) 
593 FORMAT(1X,’RA’,17,1X,'0") 
5931 FORMAT(1X, 'TX’,17,1X,F8.4) 
5932 FORMAT(1X,/TY’,17,1X,'0°) 
5933 FORMAT(1X,'TZ’,17,1X,F8.4) 
594 FORMAT(1X, ’TX’,17,1X,F8.4) 
599 FORMAT(1X,’TY’,17,1X, 0") 
600 FORMAT(1X,'TZ’,17,1X,F8.4) 
595 FORMAT(1X, ’RX’,17,1X,‘0’) 
596 FORMAT(1X,’RZ’,17,1X,'0’) 
597 FORMAT(1X,’-- BLANK LINE --’,/, 


= 1X, ‘DISPLACEMENTS APPLIED 0’) 
598 FORMAT(1X,’-- BLANK LINE --’,/, 
+ 1X, “SOLVE’,/,1X, ‘QUIT’ ) 
STOP 
END 


SUBROUTINE GAUSS (N,A) 


er. THIS SUBROUTINE PERFORMS GAUSS ELIMINATION WITH 
PIVOTING. 
INTEGER PV ! PIVOT INDEX. 
DIMENSION A(100,101) 
EPS=1.0 


10 IF (1.0+EPS.GT.1.0) THEN 
EPS=EPS/2.0 
GOTO 10 
END IF 
EPS=EPS*2.0 
EPS2=EPS*2.0 
DO 1010 I=1,N-1 
PV=I 
DOwd=14 1 
IF (ABS(A(PV,I)).LT.ABS(A(J,I))) PV=Jd 
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END DO 
IF (PV.EQ.1I) GOTO 1050 
DO JC=1,N+1 

TM=A(I,JC) 

Ei Lie) —A (PVC ) 

A(PV,JC)=TM 


END DO 
1050 IF (A(I,1I).EQ.0) GOTO 1200 ! KICK OUT IF SINGULAR 
DO JR=I+1,N ! ELIMINATION OF BELOW DIAGONAL. 


IF (A(JR,1I).NE.0.0) THEN 
R=A(JR,1I)/A(I,I) 

DO KC=I+1,N+1 
TEMP=A(JR,KC) 
A(JR,KC)=A(JR,KC) - R*A(I,KC) 

IF (ABS(A(JR,KC)).LT.EPS2*TEMP) A(JR,KC)=0.0 


C....-IF THE RESULT OF SUBTRACTION IS SMALLER THAN 2 TIMES 
& EPS TIMES THE ORIGINAL VALUE, IT IS SET TO ZERO. 

END DO 

END IF 


1060 END DO 
1010 CONTINUE 


IF (A(N,N).EQ.0) GOTO 1200 
A(N,N+1)=A(N,N+1) /A(N,N) 
DO NV=N-1,1,-1 ! BEGIN BACKWARD SUBSTITUTION. 
VA=A(NV,N+1) 
DO K=NV+1,N 
VA=VA-A(NV,K)*A(K,N+1) 
END DO 
A(NV,N+1)=VA/A(NV,NV) 
END DO | 
RETURN 


1200 PRINT*, ‘MATRIX IS SINGULAR’ 


STOP 
END 
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